Cross-diffusion induced instability on networks

The concept of Turing instability, namely that diffusion can destabilize the uniform steady state, is well known either in the context of partial differential equations (PDEs) or in networks of dynamical systems. Recently reaction-diffusion equations with cross-diffusion terms have been investigated, showing an analogous effect called cross-diffusion induced instability. In this paper, we extend this concept to networks of dynamical systems, showing that the spectrum of the graph Laplacian determines the instability appearance, as well as the spectrum of the Laplace operator in reaction-diffusion equations. We extend to network dynamics a particular network model for competing species, coming from the PDEs context. In particular, the influence of different topology structures on the cross-diffusion induced instability is highlighted, considering regular rings and lattices, and also small-world, Erd\H{o}s-R\'eyni, and Barab\'asi-Albert networks.


Introduction
A large number of real-life phenomena, for example in chemistry, ecology, and biology, give rise to a rich variety of complex behaviors, including pattern formation.The spirals that originate from chemical reactions, fish skin and animal coat patterning, and spatial vegetation patterns result from a spontaneous drive for self-organization into regular structures, both in time and space.Mathematical principles that can drive the process of pattern formation have been established by Alan Turing in 1952.In his seminal paper on the theory of morphogenesis [54], he discovered that patterns can arise as a result of the dynamical interplay between reaction and diffusion.This theory, known as Turing instability, provides a general and elegant explanation for the variety of patterns appearing in living systems: diffusion perturbs and destabilizes a homogeneous stable equilibrium, yielding to spatially inhomogeneous steady states.
On the other hand, Turing instability can also occur in networks of dynamical systems.In the seminal paper by Othmer and Scriven [41], a general mathematical framework for the analysis of instabilities in networks has been proposed and applied to regular lattices or small networks.Later it has been further explored and extended to more complex networks [4,40], and even later to multiplex [3,29,12], time-varying [45] and stochastic [5] networks.Recently, Turing bifurcations have been investigated on one-dimensional random ring networks where the probability of a connection between two nodes depends on the distance between them by using graphons [11].The surprising finding is that the conditions leading to Turing instability are the same for both, the reaction-diffusion and the network systems, but while in the continuum case, we look at the eigenvalues of the Laplace operator, in the network case we need the spectrum of the graph Laplacian.
The graph Laplacian of an undirected network with N nodes is a real, symmetric and positive semi-definite matrix, whose elements are given by where a ij are the elements of the adjacency matrix and k i is the degree of the node i, defined as a ij , a ij = 1, i ↔ j, 0 i j, i, j = 1, . . ., N.
In many areas and applications, it turns out that the eigenvalues of the graph Laplacian are useful tools, and for this reason, they have been intensively studied [2,16,27,34,37,38].The eigenvalues Λ α and eigenvectors v (α) = (v N ) ⊤ of the graph Laplacian L are determined by For an undirected network, eigenvalues are real and nonnegative.It can be easily proven that the 0 is an eigenvalue, 0 = λ 1 ≤ λ 2 ≤ • • • ≤ λ N ≤ N , and that λ 2 > 0 if the network is connected.In particular, the eigenvalue λ 2 is often called algebraic connectivity, and it holds that where |E| is the number of edges, and 2|E| can be obtained by summing the diagonal elements of the graph Laplacian.Properties of the spectrum have been investigated for particular network structures [20] including random graphs with given expected degrees [14] or more general random graphs [15].Moreover, the spectrum of the graph Laplacian, as a quantity that encodes the topological structure of the underlying graph, also influences the dynamical properties, such as the stability of the synchronous state, which has been investigated for quite some time [42,47] and has more recently been referred to as master stability function approach [9,39,44,52].The diffusion on the network given by the graph Laplacian only is a linear diffusion and considers pairwise interactions.More recently, nonlinear diffusion (also called nonlinear coupling or, in ecology, nonlinear dispersal), and higher-order interactions on networks have been proposed and investigated, see for instance [8,10,24,36].
It is widely accepted that the Turing instability is induced by diffusion and requires an activatorinhibitor scheme of interaction between agents [25].However, in the context of continuous space reaction-diffusion systems a cross-diffusion model has been proposed to describe the spatial segregation of the species [48].In this model, also known as Shigesada-Kawasaki-Teramoto (SKT) system, the reaction part does not generate an activator-inhibitor mechanism, but spatial patterns may arise thanks to cross-diffusion terms.These terms model those situations in which the movements of individuals of a species depend not only on the species itself but also on the presence of other species.This mechanism can also be observed in other contexts: for instance in the presence of chemotaxis, or also in epidemiology, where the movements of susceptible individuals are clearly influenced by the presence of infected ones [43].Then, cross-diffusion can destabilize a uniform equilibrium and induce pattern formation, even when it is not possible via standard diffusion terms; this phenomenon is known as cross-diffusion-induced instability.Reaction-cross-diffusion models have been extensively studied from different points of view and related to several applications (see [13,17,19,23,30,32,51,53] and references therein).Knowing all these ingredients, we are interested in the extension and the study of cross-diffusion and cross-diffusion-induced instability on complex networks.In addition to the natural pairing between continuous and discrete space models, this extension has been naively inspired by the discretization of reaction-cross-diffusion equations using finite differences [18].In fact, the (regular) mesh can be viewed as a regular lattice.Moreover, a natural question at this point is when/if the discretized systems show the same patterns as the continuous model (depending on the number of mesh points).
While the theory of pattern formation on networks has been developed for several network structures and dynamical rules, in this direction only a few works have been investigated [59,21], where the cross-diffusion terms considered are relatively simple (linear).We show that the theory of pattern formation on networks still holds considering general cross-diffusion terms (written in terms of the graph Laplacian).Also in this case, the conditions for cross-diffusion-driven instability are the same for the network and the continuous space case, where the eigenvalues of the graph Laplacian play the role of the eigenvalue of the Laplace operator.Then, we propose and investigate the SKT cross-diffusion model for competing species on a network.The model generalises the linear-diffusion one investigated in [49], which is known to have no heterogeneous (positive) steady states in the weak competition regime.As in the continuous setting, cross-diffusion is the key ingredient in producing non-homogeneous steady-state solutions.Our attention is however focused on the network structures (from regular rings, 2D-lattices to different complex/random graphs) in order to show, how they influence the possible dynamics of the system.In particular, we look at the spectrum of the graph Laplacian (or its distribution for random graphs).The aim of this work is mainly to point out that cross-diffusion terms can be useful ingredients in the study of complex systems, and that they can give rise to very rich dynamics depending crucially on the network topology.
The paper is organized as follows.In Section 2, we establish the general framework for crossdiffusion systems on networks and we extend the Turing instability analysis to cross-diffusion-induced instability.In Section 3 we propose a cross-diffusion model for competing species, inspired by the discretization of the SKT model on a 1D domain (presented in Appendix A).This model is analyzed on different network structures, from regular rings and lattices, to more complex networks.The non-homogeneous steady states are studied, and the distributions of the eigenvalues for different network structures are presented.The focus of this section is to understand if particular structures are more favourable to leading to cross-diffusion-induced instability.Finally, in Section 4 some concluding remarks can be found.The link between the network model and the discretization of the reaction-diffusion SKT model is made in Appendix A. The Python scripts for the simulations are freely accessible in the GitHub folder [31].

The general framework
In this section, we extend the Turing instability analysis on networks to cross-diffusion systems.This type of system can model complex natural phenomena in which nontrivial and nonlinear effects also affect the diffusion of the involved quantities.For instance, in ecology competition of species can cause migrations as a prey species tries to avoid predators, while in epidemiology susceptible individuals may try to avoid infected ones.
We consider an undirected network of N nodes.The network topology is encoded in the graph Laplacian L ∈ R N ×N .On each node, we consider two state variables (φ i , ψ i ) = (φ i (t), ψ i (t)), evolving over time.The dynamics on each node is influenced by different components: -the single node dynamics, described by the functions f, g, -standard dispersal (random movements), expressed as a diffusive flux of a species to other nodes in terms of potential difference.The diffusion coefficients are denoted by D i , i = 1, 2.
-cross-diffusion effects, due to the presence of a different species on other nodes.The crossdiffusion coefficients are denoted by D ij , i, j = 1, 2, while functions c i , i = 1, 2 describes the inter-specific competition of the two species.
-self-diffusion effects, due to the presence of the same species on other nodes.The self-diffusion coefficients are denoted by D ii , i = 1, 2, while functions s i , i = 1, 2 describes the intra-specific competition among individual of the same species.
All the considered functions are supposed to be regular enough (at least sufficiently smooth).Then the network dynamics can be modelled by where the overdot denotes time differentiation, and l ij are the elements of the graph Laplacian L as defined in equation ( 1) and they appear in the diffusion terms (standard, cross-and self-diffusion) since they encode the network structure.
As usual in Turing instability analysis, we consider a stable steady state for the single node dynamics (φ * , ψ * ), such that where the matrix J * is the Jacobian of the single node dynamics evaluated at the homogeneous steady state It turns out that the state is a homogeneous steady state for the network.Cross-diffusion induced instability arises when (φ * , ψ * ) becomes unstable to inhomogeneous perturbations.The linear stability analysis is effectively a suitable combination of the standard cases on networks and in continuous media.We introduce small perturbations δφ i , δψ i to the uniform state as and substitute this into equation (2).Linearized differential equations for δφ i , δψ i are and they can be written as where the matrix D * is the linearization of the diffusion part evaluated at the uniform steady state, given by .
We consider the spectrum of the graph Laplacian where Λ α and v (α) represent the eigenvalues and their associated eigenvectors, respectively.By expanding the perturbations δφ i , δψ i over the set of Laplacian eigenvectors v where the constants c α and b α refer to the initial conditions, system (3) is transformed into N independent linear equations for different normal modes, resulting in the following eigenvalue equation for each α: The matrix is called characteristic matrix.Then, in order to destabilize the homogeneous steady state, the characteristic matrix M * α must have a positive eigenvalue for at least one α corresponding to one eigenvalue Λ α of the graph Laplacian.Note that the eigenvalues of the graph Laplacian appear only in combination with diffusion coefficients.
Remark.Since the eigenvalues of the graph Laplacian are nonnegative, the expression of the characteristic matrix M * α turns out to be the same as in the PDEs setting, and the eigenvalues of the graph Laplacian play the role of the eigenvalues of the Laplace operator.Then the conditions on the parameters leading to the instability of the homogeneous steady state are the same.
In summary, the Turing instability analysis in the continuous setting and on discrete models can be adapted to cross-diffusion systems on networks.Note that depending on the single node dynamics, it is not always possible to destabilize the homogeneous equilibrium by means of standard diffusion terms only (in particular when the system does not have the activator-inhibitor structure).Cross-diffusion induced instability appears when the homogeneous steady state becomes unstable to inhomogeneous perturbations, due to the additional presence of cross-diffusion terms.
Remark.When Λ α = 0 the characteristic matrix (7) reduces to the Jacobian J * of the single node dynamics evaluated at the coexistence steady state.Therefore the zero-eigenvalue provides information about the dynamics of an isolated node.The non-zero eigenvalues instead account for the more complex behavior of the network.

The SKT model on networks
The goal of this section is to show that the cross-diffusion terms can destabilize the uniform equilibrium state leading to non-trivial steady states and to study the role and the influence of the network topology on these states.To achieve this, we present here a cross-diffusion model on networks for competing species, analogous to the SKT model presented in the framework of PDEs (see [13,30] and references therein).We consider two species, u and v, competing for the same resource.A system of ODEs describes the dynamics of an isolated node where r 1 , r 2 are the growth rates, a 1 , a 2 the intra-specific competition rates and b 1 , b 2 the interspecific competition rates.It is simple to show that the outcomes of ( 4) can be the total extinction, competitive exclusion or the coexistence of the species, depending on the parameter values and on the initial conditions.
We consider now a network of N nodes, its topology being fixed by the Laplacian matrix of the graph L = (l ij ) i,j=1,...,N .On each node the population sizes are denoted with u i , v i , i = 1, . . ., N and the dynamics on each node in case of isolation (neglecting the diffusion process) is given by (4).When the nodes are connected in a network, we consider diffusive movements of the individuals between linked nodes modelled in the usual way Then we consider cross-and self-diffusion effects, which cause extra movements of individuals due to intra-and inter-competition pressure.In particular, individuals living on a node leave to reach other nodes because of the presence of individuals of the same species (self-diffusion) or of the competing species (cross-diffusion).Figure 1 summarizes the gain and loss terms considered here.
In detail, they can be described by Then, the system of ODEs describing the dynamics of the network is given by where d i , D ij , i, j = 1, 2 are the coupling strengths (or diffusion coefficients).
In order to obtain cross-diffusion instability effects, it is sufficient to consider D 11 = D 22 = 0 and d 1 = d 2 = d (analogous to the PDEs case).This choice corresponds to the simpler cross-diffusion system Following [49], it can be proven that solutions with nonnegative initial conditions remain nonnegative all the time.
Proof.We must prove the fact that the set is a positively invariant region for system (6).The calculations follow [49] and it is easy to see that the extra terms coming from cross-diffusion have the correct sign configuration.
Remark.When only linear diffusion is considered [49], a priori bounds for the solutions can be obtained.With cross-diffusion terms, we cannot obtain comparison results as in [49,Theorem 3.2].
We now proceed with the linearized analysis close to the coexistence state (u * , v * ) in the weakcompetition regime (namely a 1 a 2 − b 1 b 2 > 0 and (u * , v * ) is stable for the single node dynamics (4)).Then, according to Section 2, the characteristic matrix is where Λ denotes an eigenvalue of the graph Laplacian.Remembering that the graph Laplacian has nonnegative eigenvalues, then the characteristic matrix has a negative trace and its determinant determines the stability/instability of the homogeneous steady state in which In particular, we consider the weak competition case a 1 a 2 − b 1 b 2 > 0, leading to a stable steady state where tr(J * ) < 0 and det(J * ) > 0. Note that the obtained characteristic matrix is the same as in the reaction-cross-diffusion system and the difference between the continuum and the discrete case is played by the eigenvalues of the Laplacian and the graph Laplacian respectively.The trace of the characteristic matrix is negative.Then a cross-diffusion induced instability may appear if its determinant becomes negative for at least one Λ.The determinant can be written as where the coefficients A Λ , B Λ , C Λ of the second order polynomial in d are given by In order to have a negative determinant, we need C Λ < 0. Note that there are two cases: • D 12 α + D 21 β ≤ 0, then no bifurcation can appear, independently of the network topology.
• D 12 α + D 21 β > 0, then C Λ can be negative and bifurcations can appear depending on the eigenvalues of the Laplacian of the graph.In detail, we obtain the following threshold Remark.This condition can be combined with other information on the spectrum of a particular graph to prove the stability of the homogeneous state.Since Λ N ≤ N , the simplest observation is that if N < Λ * , then no cross-diffusion induced instability is possible.We can also obtain conditions regarding the stability of a particular mode (eigenvalue).For instance, for a regular graph of N nodes and degree 2K, we know that the first nontrivial eigenvalue (also known as algebraic connectivity) satisfies Λ 2 ≤ 2KN/(N − 1).If Λ * > 2KN/(N − 1), then the mode related to Λ 2 is stable.
Remark.From equation (8) it is easy to see the cross-diffusion terms are the key ingredient to the appearance of nonhomogeneous steady states.In fact, if D 12 = D 21 = 0, then det(M Λ ).This means that, regardless of the network structure, standard diffusion terms cannot lead to Turing instability.
We look now at the determinant of the characteristic matrix as a second-order polynomial in Λ: Depending on the parameter set (the discriminant of the second-order polynomial must be positive), we can obtain an instability region Ω * = (Λ * 1 , Λ * 2 ) for the eigenvalues of the graph Laplacian.

Cross-diffusion induced instability in the SKT network
This section is mainly devoted to showing that cross-diffusion terms are able to destabilize the uniform state, leading to the so-called cross-diffusion induced instability.The key point is that this effect is not possible without them.To this end, we use a simple network topology, in order to highlight the effect per se.The set of fixed parameters appears in Table 1, while the others will be specified each time in the text.The different network structures have been generated using the Python [46] package NetworkX [26].
Remark.In the simulations, we do not vary the linear diffusion parameter d.As in the PDEs case, large values of d tend to stabilise the homogeneous state.

2K-regular ring lattice
A 2K-regular ring lattice is a graph with N nodes in a ring structure in which each node is connected to its 2K neighbors (K on either side) [57].The associated graph Laplacian is a matrix with three bands (in the centre and in the corners NE and SW), defined by The structures of the network and of the graph Laplacian are shown in Figure 2. The closed formula for the eigenvalues is 2 cos 2πk(j − 1) N , j = 1, . . ., N.
In Figure 3 the spectrum of the graph Laplacian is reported for different values of N and K, in order to study the possible appearance of patterns.In particular, if an eigenvalue is located in the instability region Ω * = (Λ * 1 , Λ * 2 ), marked with a red stripe, then system (6) admits a stable nonhomogenous steady state.We can observe in Figure 3a that the value of nonzero eigenvalues increases when K increases (and N is fixed); for small values of K the spectrum lies below the threshold value Λ * 1 while increasing K we pass from a situation in which part of the spectrum is located in the instability region to a situation in which only the smallest nonzero eigenvalue is present.Finally, for K > 25 all the nonzero eigenvalues are greater than the threshold value Λ * 2 , namely no stable pattern can appear.On the contrary, increasing N with a fixed K leads to smaller eigenvalues in the spectrum intersecting the instability region (Figure 3b).However, this will lead to different nonhomogeneous solutions.
In Figure 4 we show different outcomes of the cross-diffusion model ( 6) for the ring structure of N = 100 nodes and three different values of K, corresponding to different locations of the spectrum with respect to the instability region.Trajectories of each node over time and the final configuration of the network are reported: varying K the final configuration is not homogeneous, but a different pattern can appear.

Cross-diffusion induced instability and graph topology
This section is devoted to showing, how the cross-diffusion terms are related to the graph topology.We consider different types of graphs, including random graphs such as small world and Erdős-Rényi networks.The same parameter set describes the dynamics, but different outcomes emerge as the result of the underlying topology.As a key factor, the spectral properties of the graph Laplacian are compared.

2D-Lattices
We consider here three different two-dimensional grid graphs: triangular, square and hexagonal lattices.In a triangular lattice graph, each square unit has a diagonal edge.The square lattice has each node connected to its four nearest neighbors, while in the hexagonal lattice nodes and edges are the hexagonal tiling of the plane.These three structures constitute three different network topologies for the same set of nodes.Note that, the degree is the same for all the nodes (except for "boundary nodes"): 6 in the triangular lattice, 4 in the square lattice and 3 in the hexagonal lattice.
In Figure 5 we show the location of the spectra of the three topologies (N = 110) with respect to the instability region (obtained using the parameter set in Table 1).Note that the value of the largest eigenvalue does not change significantly increasing N , since the lattices are almost regular.It can be observed that nonhomogeneous steady states cannot appear on the hexagonal grid, while in the triangular and square lattices, the homogeneous steady state is unstable.The steady patterns on a triangular and square lattice with N = 400 nodes are shown in Figure 6 with respect to the u variable.1 and d = 0.03.Initial conditions are random perturbations of the homogeneous state.

Random graphs
We now turn our attention to several random graphs.Also in this case we want to study the influence of a particular structure on the emergence of nonhomogeneous steady states in the network.As in the previous section, we look at the nonzero eigenvalues of the graph Laplacian.Of course, since we are now dealing with random graphs, we consider the distribution of the eigenvalues (mean values and the corresponding variance) obtained with several realizations of the same random graph.
We consider the following random graphs (briefly recalling the definition and some properties).
-Regular-random graph A random-regular graph is chosen uniformly from the set of all K-regular graphs with N nodes.
-small-world graph (Watts-Strogatz) [56] Starting from a ring lattice with N nodes and K edges per node, each edge is rewired at random with probability p.This construction allows to "tune" the graph between K-regular rings (p=0) and random graphs (p = 1).The small-world region (10 −4 < p < 10 −1 ) is characterized by a small average path length and a large clustering coefficient.
-binomial graph (Erdős-Réyni ) [22] Considering N nodes, each of the possible edges is chosen with probability p.
-preferential attachment model (Barabási-Albert) [6] A graph of N nodes is grown by attaching new nodes each with K edges that are preferentially attached to existing nodes with a high degree.
We compare these different graph topologies and the appearance of cross-diffusion-induced instability.As in the previous sections, we are interested in the spectrum of the graph Laplacian.Since we deal with particular classes of random graphs, general results cannot be achieved by just looking at a particular realization.For each type of structure, we fix the number of nodes N = 100 varying the other network parameters (the probability p for Watts-Strogatz and Erdős-Réyni graphs, the number of edges of a node for random-regular and Barabási-Albert graphs).We generate 1000 realizations for each structure and each parameter value, obtaining a distribution of the eigenvalues.In Figure 8, the mean (dots) and the variance (error bars) of the eigenvalues of the graph Laplacian are shown for the four types of structures.The first nontrivial eigenvalue and the last one are marked in green.The red stripe denotes the instability region (Λ * 1 , Λ * 2 ), while the green solid horizontal line denote the threshold Λ * , relevant to the parameter set in Table 1.
In a regular-random graph (Figure 7a) it can be seen that the cross-diffusion induced instability appears on average in the region when 2 ≤ K/2 ≤ 13 and N = 100, while greater values of the nodedegree only lead to a very small probability of instability and the uniform steady state is expected to be stable.In a small-world network, with K/2 = 15 being fixed, we report in Figure 7b the averages and variances of the eigenvalues for different values of the probability p of rewiring an edge in a larger interval than the small-world regime.It turns out that the systems show with a probability well bounded away from zero non-uniform steady states for the considered values of p because at least one averaged eigenvalue is always located well within the instability region.Depending on the number of eigenvalues in the instability region, different types of non-uniform steady states can be observed.For the Erdős-Rényi graph, cross-diffusion induced instability is most likely to appear for small values of the probability p of choosing an edge, while it is not likely to occur for larger values.In order to better compare this structure with the other types, we consider p ≈ 4K/(N − 1), which for large N approximates the average number of links for one node.With N = 100 and K = 30, we obtain p ≈ 0.3, for which the systems display cross-diffusion induced instability with a probability well bounded away from zero.Finally, the Barabási-Albert graph is quite different from the others.The A particular realization of the random graph (N = 100) is reported, in which the edges between "neighbors node" (K/2 = 15 on each side) appear in magenta, while "long-range interaction" in green.
spectrum has a large variance, and especially for large values of K (that in this case represents the number of edges chosen for every new node in the growing process).We observe that the possibility of pattern formation really depends much more on the particular realization of the graph in comparison to the other graph structures.
A final observation has to be made on the non-uniform steady state that arises through crossdiffusion induced instability.Its shape is still quite regular and "smooth" as reported for regular rings in the previous section, but small variations from this regular configuration appear (see Figure 9).For the other structures instead, the system tends to a stable non-uniform configuration that however does not present the characteristic shape of valleys and bumps.
In addition to the phenomenological discussion, ecological back-interpretation is also important.The first crucial observation is that cross-diffusion enables the coexistence of the species with a nonhomogeneous distribution.In fact, without cross-diffusion, the homogeneous distribution of the populations (namely all nodes reach the same steady state) is the only outcome.This fact is indeed important in ecology: homogeneity and synchronization of metacommunities may be harmful to species' survival.Cross-diffusion also affects the populations' abundance as observed in the PDEs model.Looking at the total abundance of the species on the network for the simulations presented in the papers and the parameter set considered, we observe a small decrease in the total abundance of population u and an increase in the total abundance of population v with respect to the homogenous case (without cross-diffusion).For instance, for the ring topology (N = 100, K = 10), parameter set as in Table 1 and d = 0.03, species u decreases of 1.78% while species v increases of 12.3%.This quantitative difference reflects the difference in the order of magnitude of the population sizes at the homogeneous steady state (see for instance Figure 4) and it is a limitation of the particular parameter set.However, the qualitative trend highlights a counter-intuitive effect driven by cross-diffusion, namely that although species u tries to avoid v, it is actually v that benefits from this.However, it is worthwhile to say that the total population abundance is not the only factor measuring the advantages/disadvantages of the two competing populations.), while the green solid horizontal line denotes the threshold Λ * , relevant to the parameter set in Table 1 and d = 0.03.

Conclusion and Outlook
In this paper, we have extended the cross-diffusion induced instability, already known in reactioncross-diffusion systems of PDEs, to networks of dynamical systems.The non-standard diffusion terms are also written using the graph Laplacian.We established the general framework and we have shown that in this case the linearization slightly differs from the standard case and the PDEs case.Then, we adapted and investigated the SKT cross-diffusion model for competing species on a network.As already known in the context of reaction-cross-diffusion models, the cross-diffusion terms can destabilize the homogeneous equilibrium state and cause the appearance of organized states and patterns, which is  not possible only with standard diffusion terms.The important finding is that the obtained conditions for cross-diffusion induced instability in the network framework are the same as the continuous case, where the role of the eigenvalues of the Laplace operator is played by the eigenvalues of the graph Laplacian.Also in this case, an instability region depending on the model parameter characterizes the possibility to observe patterns, and the location of the eigenvalues of the graph Laplacian determines their presence and shape.Finally, we have analyzed different network structures (such as regular rings, 2D lattices and different random graphs) in order to show how the network structure influences the possible outcomes of the system.In particular, we have looked at the spectrum of the graph Laplacian (or its distribution for random graphs).We conclude that cross-diffusion induced instability depends on the network structure: for instance, for the SKT network model in the weak competition case it is more likely to appear on triangular lattices, and on small-world and Barabasi-Albert random graphs.
One key aim of this work is to point out that cross-diffusion terms can be a useful tool to model complex systems, and that they can give rise to richer dynamics.Several research directions arise at this point.On the one hand, we can look at the dynamical and bifurcation aspects of this topic.As widely investigated for cross-diffusion PDEs systems, we want to find an entropy functional or a Lyapunov function for networks in order to achieve global stability results [7,35,50].Furthermore, a deeper investigation of the system outcomes combined with the bifurcation structure (that can be computed by the continuation software pde2path [55]) may reveal the presence of periodic patterns, as in the PDEs case.On the other hand, one could take an even more detailed look at the influence of the graph topology.In this regard, the parallel between continuous and discrete models is intriguing.Inspired by the discretized version of the PDEs model involving the Laplace operator, we considered in this paper different graph topologies, changing completely the corresponding operator in the continuous case.Therefore, it would be interesting to understand the limit model back to the continuous case.Moreover, following [4], it seems possible to extend the analysis to directed networks.Moreover, considering for instance the SKT model on networks, it can also be possible to consider the dynamics of the two competing species evolving on different networks, as sketched in Figure 10.This would lead us to consider the dynamics on a multilayer network [1,12,28], where the two layers can be different, for which the extension of the theory of Turing patterns has already been presented in [3,12,29].It could also be interesting to study other types of cross-diffusion terms involving different nodes: for instance, in the SKT model, the movement of one species from node i to node j is influenced by the presence of the same species or the competing one on node j.In this scenario, there is an information flow or knowledge about the status of the other nodes, or it can be seen as a weighted network with link weights depending on the quantities on each node.Finally, one could study possible applications of cross-diffusion systems on networks, ranging from ecology and landscape modelling to disease spreading [21,33], that are characterized by a strong interplay between dynamics and network structure.

A Discretized SKT reaction-diffusion model
The extension of cross-diffusion systems on networks comes from the discretization of reaction-crossdiffusion systems by finite differences.In particular, we consider the SKT model with self-and cross-diffusion terms given by where the quantities u(t, x), v(t, x) ≥ 0 represent the population densities of two species at time t and position x, confined and competing for resources on a bounded and connected domain Ω ⊂ R n .The parameters are the same as the network model presented in Section 3.
The discretized SKT model on a 1D domain of length ℓ using N nodes and uniform mesh size h = ℓ/(N − 1) is [18] The uniform mesh of N nodes can be seen as a 1D lattice of N nodes, each one connected to its two nearest neighbors (where the boundary conditions determine the edges of the boundary nodes).Indicating with l ij , i, j = 1, . . ., N are the elements of the graph Laplacian, the discretized system can be written in the form l ij u j v j , i = 1, . . ., N. (11) Note that the graph Laplacian L ∈ R N ×N is the usual tridiagonal matrix with eigenvalues given by the closed formula [58] Λ j = 2 − 2 cos π(j − 1) N , j = 1, . . ., N.
The main difference with the network model ( 5) is the factor 1/h 2 in front of the diffusive terms coming from the finite difference scheme.Roughly speaking, it allows the eigenvalues of the graph Laplacian to approximate the first N eigenvalues of the Laplace operator when h becomes small (namely increasing the number of mesh points/nodes).

Figure 1 :
Figure 1: Schematic representation of the dynamics on a node in the network.

r 1 r 2 a 1 a 2 b 1 b 2 d d 12 d 21 Table 1 :
Set of parameter values of the single node dynamics and diffusion coefficients of the SKT model (6) used in the numerical simulations.The set r i , a i , b i , (i = 1, 2) corresponds to the weak competition case (a 1 a 2 − b 1 b 2 > 0), namely the homogeneous steady states is stable for the single node dynamics.The remaining parameters related to the network structure will be specified in the text.

Figure 2 :
Figure 2: Structure of a 2K-regular ring lattice and of its graph Laplacian.(a) Structure of the regular 2K-ring with N = 30, K = 8, in which the 2K nearest neighbors (K on each side) of a particular node are highlighted in magenta.(b) Structure of the corresponding graph Laplacian, where grey dots indicate non-zero elements.

Figure 3 :
Figure 3: Cross-diffusion driven instability on regular rings.Grey dots mark the eigenvalues, while the first and the last nonzero eigenvalues appear in green.The instability region Ω * = (Λ * 1 , Λ * 2 ) and the threshold Λ * are marked with the red stripe and the green horizontal solid line, respectively (relevant to the parameter set in Table 1 and d = 0.03).(a) Spectrum of the graph Laplacian of the regular 2K-ring with N = 100 and different values of K. (b) Spectrum of the graph Laplacian of the regular 2K-ring with K = 10 and different values of N .

Figure 4 :
Figure 4: Dynamics on a 2K-regular ring of N = 100 nodes for different values of K.The parameter set is reported in Table 1 and d = 0.03.Initial conditions are random perturbations of the homogeneous state.Purple and orange correspond to u and v, respectively.(Left) Dynamics over time.The trajectory of one node is highlighted, while the others are shown in a lighter color; each node is identified via its index i ∈ {1, 2, . . ., N }.(Right) Stable configurations on the ring.Different kinds of patterns appear.11

Figure 5 :
Figure 5: Cross-diffusion driven instability on triangular (△), square ( ) and hexagonal ( ) lattices (with N = 110).Grey dots mark the eigenvalues, while the first and the last nonzero eigenvalues appear in green.The instability region Ω * = (Λ * 1 , Λ * 2 ) and the threshold Λ * are marked with the red stripe and the green horizontal solid line, respectively (relevant to the parameter set in Table 1 and d = 0.03).

Figure 6 :
Figure 6: Steady state (u variable) on a triangular and square lattices with N = 400 nodes.The parameter set is reported inTable 1 and d = 0.03.Initial conditions are random perturbations of the homogeneous state.

Figure 7 :
Figure7: Comparison of different graph topologies (regular-random, small-world, Erdős-Rényi and Barabasi-Albert).A particular realization of the random graph (N = 100) is reported, in which the edges between "neighbors node" (K/2 = 15 on each side) appear in magenta, while "long-range interaction" in green.

Figure 8 :
Figure8: Comparison of different graph topologies and the appearance of cross-diffusion induced instability.The mean and variance of the eigenvalues of the graph Laplacian, obtained with 1000 realization of the random graphs are denoted by dots and error bars, respectively; the first nontrivial eigenvalue and the last one are marked in green.The red stripe denotes the instability region (Λ * 1 , Λ * 2 ), while the green solid horizontal line denotes the threshold Λ * , relevant to the parameter set in Table1and d = 0.03.

Figure 9 :
Figure 9: Non-homogeneous steady state on a small-world network with N = 100 and K = 30 for different values of the probability p.The parameter set is reported inTable 1 and d = 0.03.Initial conditions are random perturbations of the homogeneous state.

Figure 10 :
Figure 10: Different way to model the dynamics of two species on a network: the species evolve, compete and move (a) on the same network, or (b) on a multilayer network.

Table 1
and d = 0.03.Initial conditions are random perturbations of the homogeneous state.Purple and orange correspond to u and v, respectively.(Left) Dynamics over time.The trajectory of one node is highlighted, while the others are shown in a lighter color; each node is identified via its index i ∈ {1, 2, . . ., N }.(Right) Stable configurations on the ring.Different kinds of patterns appear.11 Table 1 and d = 0.03.Initial conditions are random perturbations of the homogeneous state.