The linear noise approximation for reaction-diffusion systems on networks
Stochastic reaction-diffusion models can be analytically studied on complex networks using the linear noise approximation. This is illustrated through the use of a specific stochastic model, which displays traveling waves in its deterministic limit. The role of stochastic fluctuations is investigated and shown to drive the emergence of stochastic waves beyond the region of the instability predicted from the deterministic theory. Simulations are performed to test the theoretical results and are analyzed via a generalized Fourier transform algorithm. This transform is defined using the eigenvectors of the discrete Laplacian defined on the network. A peak in the numerical power spectrum of the fluctuations is observed in quantitative agreement with the theoretical predictions.
pacs:89.75.Kd, 89.75.Fb, 05.10.Gg, 02.50.-r
Pattern formation in reaction-diffusion models is a subject which finds applications in many distinct fields, including ecology mimura (); maron (); baurmann (); riet (), biology meinhardt (); harris (); maini (); bhat (); miura () and chemistry prigogine (); castets (); quyang (). In his seminal paper, Turing turing () pointed out that a small perturbation of a homogeneous state in a reaction-diffusion system can, due to an instability, rapidly grow to eventually yield stable non-homogeneous patterns. The vast majority of studies devoted to investigating the emergence of Turing patterns use deterministic partial differential equations to model the reactions and diffusion of the constituents, which are characterized by continuous spatial distributions.
These ideas have recently been generalized in two distinct directions. First, the concept of a Turing instability can be extended to embrace systems defined on complex networks nakao (). This is an important step forward vespignani (), which should eventually shed novel light onto the mechanisms that drive self-organization on networks colizza (); satorras (); colizzavespignani (). Second, Turing patterns have also been observed and analyzed in models with a finite number of constituents. In this case the reaction-diffusion processes are described by individual-based models which take into account the intrinsic discreteness of the system. Stochastic effects are therefore present, and ultimately stem from the finite size of the population of elementary constituents. In this paper we build on earlier work ourwork (), to bring these two features together, and study stochastic patterns on networks. Specifically, we predict the existence of stochastic waves on networks, and show how they may be described analytically. We illustrate this on a specific model, but it should be clear from the analysis that we present that the ideas and techniques are generally applicable.
At first sight it might appear surprising that stochastic effects are important in reaction-diffusion systems, which after all consist of a large number of constituents. However, the fluctuations due to the discreteness of these constituents can amplify through resonant effects and so yield macroscopically ordered patterns, both in time mckanePRL (); dipatti () and in space goldenfeld (); biancalani (); woolley (). Stochastic Turing patterns biancalani () (also termed quasi-Turing patterns goldenfeld ()) can appear in a region of the parameter space for which the homogeneous fixed point is predicted to be stable from a deterministic linear stability analysis. Similarly, stochastic waves tommasoOnde () have been observed in reaction-diffusion models defined on a regular lattice. Importantly, the effect of fluctuations arising from this discreteness can not only be seen in numerical simulations, but also analytically understood by expanding the governing master equation within the so-called Linear Noise Approximation (LNA) scheme.
Starting from this setting, in ourwork () stochastic Turing patterns have been predicted, and numerically observed, on a scale free network, so extending the conclusions of nakao (), beyond the deterministic setting. Here we report analytical progress by carrying out the LNA for systems in which the population is distributed over a set of nodes, which are connected to each other in some way. Each node hosts a large number of individuals, so this is a metapopulation model hanski () — a collection of populations allowing some exchange between them. Examples are: in ecology, where individuals reside on patches and may migrate to other patches that are nearby hanski (); in island models of evolutionary theory, where individuals carrying certain alleles may migrate to other islands maruyama (); in epidemiology, where the nodes are cities connected by commuters who carry disease ganna () and in reaction kinetics, where the nodes are compartments in which chemical reactions take place joe ().
More concretely, we shall consider a stochastic version of the Zhabotinsky model zhab (), introduced into a static scale-free network. This latter is created via the preferential attachment probability rule barabasi (). The power spectrum of fluctuations is analytically calculated by developing and systematizing the LNA technique to network-based applications. A localized peak for the power spectrum signals the presence of stochastic travelling waves, a prediction that we confirm with stochastic simulations gillespie (). The power spectrum is calculated from a generalized Fourier transform, the standard plane waves found in a spatial context being replaced by the eigenvectors of the discrete Laplacian operator defined on the network. To benchmark theory and simulations, we have therefore implemented and tested a numerical routine which handles the generalized Fourier analysis. This is a diagnostic tool that could prove useful beyond the specific case study, by guiding the unbiased search for structured patterns on a network topology kouvaris ().
The paper is organized as follows. In the next section we will introduce the stochastic model and discuss the basic steps of the LNA analysis on a graph. The technical aspects of the presentation will be relegated to specific Appendices. In Section III the linear stability analysis for the model in its deterministic limit is carried out. Numerical simulations are performed to show that the deterministic wave manifests itself as a localized peak in the power spectrum, as obtained from the generalized Fourier transform. In Section IV we derive the power spectrum of fluctuations which points to the existence of traveling waves, seeded by inherent stochasticity, outside the region of deterministic instability. Stochastic simulations confirm the validity of the theory. In the final Section we sum up and conclude.
Ii Model definition and the linear noise approximation (LNA)
The reaction scheme that we will investigate was introduced by Zhabotinsky et al. to study travelling waves arising from destabilization of a homogeneous state zhab (). The scheme involves molecules of three chemical species: , and — that we will also respectively call the first, second and third species. The molecules are placed on the nodes of a network composed of nodes, each of which has a finite volume . We label a molecule of species located on the -th node with ; and are similarly defined. The number of molecules of type , and are denoted by , and , respectively. The -dimensional vectors: , and , specify the state of the system. Within each node, the molecules interact through the following reaction scheme:
The reaction rates are denoted by . As explained in zhab () they are all constant except that is given by , with .
The structure of the network is described by the adjacency matrix, . This is a symmetric matrix whose elements, , is equal to one if node is connected to node , and zero otherwise. The molecules can migrate between two connected nodes as specified by the diffusion reactions:
The constants , and are the diffusion coefficients.
The construction of a stochastic model proceeds by assigning a transition rate to each reaction. They indicate the probability per unit of time to transit from state to state . To lighten the notation, we only write the components of the vectors which refer to molecules that take part in a reaction in the transition rates. Invoking mass action, the transition rates associated with reactions (II) read vanKampen ():
In a similar way the transition rates for the diffusion reactions (2) are given by:
As the dynamics is Markovian, the probability density function (PDF) that the system is in state at time , , satisfies the master equation:
This is the fundamental equation that governs the dynamics of the system.
The LNA can be applied by carrying out the van Kampen expansion for the master equation vanKampen (). This begins with changing variables from to , where :
The functions , and describe the concentrations of each chemical species in the deterministic limit, that is, obtained by letting . In this limit the system is not subject to fluctuations. As shown in Appendix A, the deterministic concentrations satisfy the following system of ordinary differential equations:
Hereafter, a dot above a symbol indicates the time derivative taken with respect to the rescaled time, . The symbol denotes the Laplacian operator defined on the network and reads:
where is the connectivity of node , . This form of the Laplacian operator vulpiani () reflects our specific choice for the microscopic reaction rates (II). Other choices are possible which would yield modified Laplacians sneppen (). Working in the context of the proposed formulation, is symmetric, a feature of which we shall take advantage of when performing the generalized Fourier analysis described below.
For finite volume , the system is subject to intrinsic noise; these fluctuations perturb the solutions of the deterministic model (II), which describes the dynamics of the model in the limit . Within the LNA, the fluctuations are Gaussian and given by a linear Fokker-Planck equation. Before turning to discuss the role played by stochastic fluctuations, we will start by focusing on the deterministic scenario. We will in particular derive the conditions under which self-organized patterns of the wave type emerge. The next section is devoted to this issue.
Iii Pattern formation in the deterministic limit
The analysis of pattern formation for system (II) defined on a regular lattice in the continuum limit has been already carried out in zhab (). Here, we review some of the results of zhab (), before moving on to discuss how the network affects the pattern formation. Throughout our analysis we have used a scale-free network generated with the Barabási-Albert preferential attachment algorithm barabasi (), with nodes and mean degree .
We first establish contact with the notation adopted in zhab () by making the following choices: , , , and . As in zhab (), we fix some of the parameters: , and . We also set . The parameters can be freely adjusted and select different dynamical regimes.
The system of differential equations (II) admits three fixed points zhab (). One of these corresponds to the extinction of both and species and is always stable. Another one is a saddle. The last one is non-trivial, and its stability depends on the values of . It is around this point in the two-dimensional plane defined by and that the pattern formation is investigated. The concentrations at the fixed point are independent of () and can be numerically determined.
Patterns arise when becomes unstable with respect to inhomogeneous perturbations cross (). To look for instabilities, we introduce small deviations from the fixed point, , and linearise system (II) around it:
The explicit expressions for the matrices and are given in Appendix A (the label NS stands for “non-spatial” and SP for “spatial”). For a regular lattice, the Fourier transform is usually employed to solve the above linear equations. This analysis needs to be adapted in the case of a system defined on a network. To this end we follow the approach of nakao (); ourwork () and start by defining the eigenvalues and eigenvectors of the matrix :
Since the Laplacian is symmetric, the eigenvalues are real and the eigenvectors form an orthonormal basis. It can actually be proven that for a case of a Barabási-Albert network the are negative and non-degenerate nakao (); ourwork (). We can now define a transform based on the eigenvectors which takes the role that the Fourier transform took on for a regular lattice. This leads to the following transforms which will be used throughout the remainder of the paper:
where is any function of the nodes and of time. This is a standard Fourier transform in time, but with the spatial Fourier modes replaced by the eigenvectors of the network Laplacian. If the network is a regular lattice, the transform (III) reduces to a standard Fourier transform for discrete space. From now on the index is used to label the variable conjugate to the nodes. Using this definition, one can define the power spectrum as .
that is now decoupled in the nodes and in time and thus readily solvable. The matrix for a given , is a matrix whose eigenvalues characterize the response of system (II) to external perturbations. The eigenvalue with the largest real part will be denoted by . If the fixed point is unstable and the system exhibits a pattern whose spatial properties are encoded by . This is the analog of the wavelength for a spatial pattern in a a system defined on a regular lattice; it is customarily written in this case. When the imaginary part of the eigenvalue, , is different from zero, the pattern oscillates in time cross (). A system unstable for and is said to undergo a wave-instability and the emerging patterns consist of travelling waves.
In Fig. 1, left panel, the domain of instability is shown as a shaded region in the plane . The fixed point is stable for fixed when . At a wave instability sets it and travelling waves are found to occur for . The real and imaginary parts of the eigenvalues are depicted in the right panel, as a function of , for two choices of the parameters , for which the system is respectively stable and unstable. The circles in the left panel of Fig. 1 indicate these two choices.
Since the system is defined on a network, the emerging patterns present two main differences as compared to those obtained for the case of conventional reaction-diffusion models defined on the continuum. First, only some of the wavelengths, , are allowed. This is due to the fact that the solutions of Eq. (10) form a discrete set; such a feature also occurs for systems defined on periodic lattices. In this latter case however, the wavelengths are equally spaced and proportional to the lattice spacing. By contrast, for systems defined on a complex network, there is no clear periodic structure and the wavelengths are clustered or irregularly distributed, as displayed in panel (b) of Fig. 1. The second unusual trait has to do with the shape of the patterns. In a reaction-diffusion system defined on a regular lattice, each point of of the lattice has a concentration which assumes a value significantly different from that characterising the homogeneous fixed point. However, for a network, only a fraction of nodes have concentrations which are significantly different from that of the homogeneous fixed point. The fraction that are differentiated in this way depends on the connectivity of the nodes and on the ratio of the diffusion coefficients nakao (). This feature cannot be simply understood from linear stability analysis, as it relates to the localization of the Laplacian eigenvectors in large networks, a property that has been recently investigated in this context in menzinger ().
In Fig. 2 the power spectrum of the concentration is plotted for a choice of the parameters that correspond to the leftmost circle (blue online) in Fig. 1(a). A peak is displayed for () (), in complete agreement with the predictions of the linear stability analysis. Similar results are obtained for the other concentrations and . Thus, the generalized Fourier algorithm based on Eqs.(III) can be effectively employed to resolve complex patterns that develop on networks. This is a valuable tool which, we believe, could prove useful for the many applications where the dynamics on a network is well known to be central, from neuroscience to epidemics.
The next section is entirely dedicated to the study of stochastic patterns, aiming at generalizing the deterministic picture of Fig. 1. By applying the LNA, we will demonstrate that stochastic waves exist in a region of the parameter space for which the deterministic analysis predicts a stable homogeneous fixed point. The presence of stochastically driven patterns will be revealed by an analytical calculation of the power spectra of fluctuations. The theoretical predictions will then be validated by reconstructing the power spectrum from the stochastic time series. Properties of the patterns that in the deterministic picture depend on the eigenvectors, such as localization, will not be addressed in the present work.
Iv Power spectra of fluctuations and stochastic patterns
While in the deterministic limit a study of the eigenvalues reveals the range of parameter values for which patterns are expected to occur, this prediction is not conclusive for systems that are subject to noise. Simulations of the master equation have shown that patterns arise even for parameter values for which the underlying fixed point is stable, provided that the system is sufficiently close to an instability. The corresponding patterns have been called stochastic patterns biancalani () or quasi-patterns goldenfeld (). The LNA allows one to gain analytical insight into the mechanism that yields such patterns lugo (); deanna (); biancalani (); goldenfeld (). The method has been extended in ourwork () to the case of a reaction-diffusion system on a network. Here we shall develop the method further and provide the first evidence for the spontaneous emergence of stochastic waves (or quasi-waves) on a network.
The stochastic perturbations about the fixed point are described by the variables ( and ) as defined in (LABEL:eq:ansatz) with , and .
As discussed in Appendix A, the role of fluctuations can be quantified by use of the van Kampen system-size expansion, which is equivalent to assuming the LNA. The quantity acts as the small parameter in the perturbative expansion: at the leading order, the macroscopic deterministic equations (II) are recovered. At the next-to-leading order one obtains the Fokker-Planck equation (18) for the distribution of stochastic fluctuations . To quantify the impact of the stochastic components of the dynamics, the power spectrum of fluctuations is utilized. It is defined as:
where is the complex conjugate of , which is given by transforming via the generalized transform (III). The average is performed over many realizations of the stochastic dynamics. Following the strategy outlined in Appendix B, one can then derive an analytic expression for the power spectrum of fluctuations. This is equation (B). Once the parameters of the model have been assigned, it is therefore possible to calculate the power spectrum of fluctuations and look for signatures of emerging self-organized structures. In Fig. 3 the analytical power spectrum for species is plotted for a choice of parameters that corresponds to the rightmost circle in Fig. 1 (a), namely outside the region for which the deterministic waves occur.
As can be seen, the power spectrum of fluctuations is characterized by a localized peak for . Therefore, species oscillates with an angular frequency and, at the same time, displays a pattern at wavelength . Stochastic waves, or quasi waves, are hence predicted to occur, in a region of the parameter plane for which the homogeneous fixed point is stable, according to the deterministic picture. In other words, stochastic corrections, stemming from finite size, and, as such, endogenous to the system under scrutiny, can eventually produce macroscopically ordered structures.
To test the correctness of the theoretical prediction we carried out stochastic simulations of the processes (II) and (2) using the Gillespie algorithm gillespie (). The numerical power spectrum is reconstructed by applying the generalized transform (III) to the time series, and averaging over independent realizations of the stochastic dynamics. The result is shown in Fig. 4 and is seen to agree with the theoretically-predicted spectrum. The location of the maximum is captured by the theory, as well as the characteristic shape of the profile.
Pattern formation has been extensively studied in the literature and with reference to a wide variety of problems. Typically, the concentrations of the species involved are assumed to obey partial differential equations. The conditions under which an instability occurs follow from standard linear stability analysis around a stable homogeneous fixed point. Recently, the emergence of steady state inhomogeneous patterns has been also studied for deterministic reaction-diffusion models defined on a network, generalizing the concept of a Turing instability to this important new area.
Deterministic models represent, however, an idealized approach to the phenomenon being investigated: they omit stochastic fluctuations that need to be included when dealing with finite populations of interacting elements. Finite size corrections result in intrinsic stochastic perturbations which undergo amplification, and through a resonant mechanism eventually yield self-organized patterns.
In this paper, we have further developed the theory of stochastic patterns for reaction-diffusion systems defined on a network. The analysis is based on a systematic and general application of the linear noise approximation scheme. Stochastic traveling waves are predicted to exist, and numerically observed in a region of the parameter plane for which the deterministic set of partial differential equations converges to a stable homogeneous solution. The analysis is carried out for a specific system, the stochastic analog of the deterministic model introduced in zhab (). The techniques discussed are however general, and can be readily adapted to any reaction diffusion model defined on a network. To benchmark theory and simulations we have also developed, and successfully tested, a numerical algorithm that performs the generalized Fourier transform, employed in the analytical derivation. This transform decomposes the signal along the eigenvectors of the discrete Laplacian operator, tailoring the analysis to the network under consideration, so allowing the spectral properties of the emerging patterns to be fully characterized.
Acknowledgements.T. B. thanks the EPSRC (UK) for partial financial support. D.F. acknowledges financial support from Ente Cassa di Risparmio di Firenze and the Program Prin2009 funded by MIUR (Italy).
Appendix A Van Kampen system-size expansion
In the main text, the van Kampen system-size expansion was used to approximate the master equation (5) by a deterministic system of ordinary differential equations — that describes the macroscopic evolution of the concentrations — together with a linear Fokker–Planck equation which characterizes the fluctuations around the macroscopic solution. This Appendix details the calculations using the system-size expansion.
In the following (and throughout the paper) the indexes and refer to the nodes of the network and range from to . The indexes and label the chemical species and range from one to three. Finally, the -dimensional vectors, such as , and , are displayed in bold.
Master equations, such as (5), can be rewritten by making use of step operators, , which represent the creation/destruction of a molecule of species at node . For instance, for species , they act on a general function by
The master equation (5) then reads:
We now apply the change of variable (LABEL:eq:ansatz). In the new variables, the step operators admit an expansion in powers of vanKampen (). The LNA corresponds to the truncation of the expansion at second order, namely:
The left-hand side of Eq. (A) can be expressed in terms of the PDF of the new variables, . This implies that
In the above equation there are terms which are either or . By contrast, the right-hand side of Eq. (A) contains and terms. They can be balanced by rescaling time by . Collecting together the terms of the same order and setting their sum at each order to zero gives, at the leading order, the deterministic system (II). Likewise, the next-to-leading order yields the Fokker-Planck equation for the fluctuations:
The form of the matrices and follow from the expansion of the transition rates (II) and (II). For illustrative purposes, here we shall discuss only the first of the transition rates (II), , explicitly. The contribution to matrix associated with this term, labeled , is found to be
Clearly, the only non-zero entry is for , since the rate involves only the species. The other diffusion rates, and , yield similar contributions for respectively and , with diffusion coefficients and concentrations corresponding to the diffusing species. The contributions arising from the transition rates for the reactions (II) follow in a similar fashion.
In most applications, the main point of interest is to study the fluctuations around a fixed point. This is certainly so our case, as we aim to characterize the pattern that originates from a small perturbation of the fixed point . We therefore substitute , and and label by and the matrices evaluated at the fixed point.
From the form of the reaction rates it is clear that the following decompositions hold ourwork ():
We end by giving the elements of the matrix and . The elements of are
and those of matrix are
Appendix B Analysis of the Power Spectra
The Gaussian white noises have zero mean and correlator:
Equation (B) generalises (9) to include stochastic fluctuations. In solving Eq. (B), we again make use of the transforms (III). We express the and the associated noise in terms of their transformed analogs. Collecting each term, except the noise, to the left-hand side of the equation yields:
where is the identity matrix. By introducing the solution of Eq. (25) may be written as:
The symbol signifies the adjoint operator, here equivalent to the conjugate transpose operator. We now need to express in terms of known quantities. We begin by transforming Eq. (24) using the inverse transform (III), which leads to
The dependence on the Laplacian eigenvectors can be eliminated using the fact that they are orthonormal and complete:
- (1) M. Mimura and J. D. Murray, J. Theor. Biol. 75, 249 (1978).
- (2) J. L. Maron and S. Harrison, Science 278, 1619 (1997).
- (3) M. Baurmann, T. Gross and U. Feudel, J. Theor. Biol. 245, 220 (2007).
- (4) M. Rietkerk and J. van de Koppel, Trends Ecol. Evol. 23, 169 (2008).
- (5) H. Meinhardt and A. Gierer, BioEssays 22, 753 (2000).
- (6) M. P. Harris, S. Williamson, J. F. Fallon, H. Meinhardt and R. O. Prum, Proc. Natl Acad. Sci. USA 102, 11734 (2005).
- (7) P. K. Maini, R. E. Baker and C-M. Chuong, Science 314, 1397 (2006).
- (8) S. A. Newman and R. Bhat, Birth Defects Res. (Part C) 81, 305 (2007).
- (9) T. Miura and K. Shiota, Dev. Dyn. 217, 241 (2000).
- (10) I. Prigogine and R. Lefever, J. Chem. Phys. 48, 1695 (1968).
- (11) V. Castets, E. Dulos, J. Boissonade, and P. De Kepper, Phys. Rev. Lett. 64, 2953 (1990).
- (12) Q. Ouyang and H. L. Swinney, Nature 352, 610 (1991).
- (13) A. M. Turing, Phil. Trans. R. Soc. Lond. B 237, 37 (1952).
- (14) H. Nakao and A. S. Mikhailov, Nature Physics 6, 544 (2010).
- (15) R. Pastor-Satorras and A. Vespignani, Nature Physics 6, 480 (2010).
- (16) V. Colizza, A. Barrat, M. Barthelemy and A. Vespignani, Proc. Natl Acad. Sci. USA 103, 2015 (2006).
- (17) R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001).
- (18) V. Colizza, R. Pastor-Satorras and A. Vespignani, Nature Physics 3, 276 (2007).
- (19) A. J. McKane and T. J. Newman, Phys. Rev. Lett 94, 218102 (2005).
- (20) T. Dauxois, F. Di Patti, D. Fanelli and A. J. McKane, Phys. Rev. E 79, 036112 (2009).
- (21) T. Butler and N. Goldenfeld, Phys. Rev. E 80, 030902(R) (2009).
- (22) T. Biancalani, D. Fanelli and F. Di Patti, Phys. Rev. E 81, 046215 (2010).
- (23) T. E. Woolley, R. E. Baker, E. A. Gaffney and P. K. Maini, Phys. Rev. E 84, 046216 (2011).
- (24) T. Biancalani, T. Galla and A. J. McKane, Phys. Rev E 84, 026201 (2011).
- (25) M. Asslani, F. Di Patti and D. Fanelli, Phys. Rev. E 86, 046105 (2012).
- (26) I. Hanski, Metapopulation Ecology (Oxford University Press, Oxford, 1999).
- (27) T. Maruyama, Stochastic Problems in Population Genetics. Lectures in Biomathematics 17 (Springer, Berlin, 1977).
- (28) G. Rozhnova, A. Nunes and A. J. McKane, Phys. Rev. E 84, 051919 (2011).
- (29) J. D. Challenger and A. J. McKane. arXiv:1302.0362.
- (30) A. M. Zhabotinsky, M. Dolnik and I. R. Epstein, J. Chem. Phys. 103, 10306 (1995).
- (31) A-L. Barabási and R. Albert, Science 286, 509 (1999).
- (32) D. T. Gillespie, J. Comp. Phys. 22, 403 (1976).
- (33) N. E. Kouvaris, H. Kori and A. S. Mikhailov, PLoS ONE 7, e45029 (2012).
- (34) N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North Holland, Amsterdam, 1992).
- (35) R. Burioni, S. Chibbaro, D. Vergni and A. Vulpiani, Phys. Rev. E 86, 055101(R) (2012).
- (36) I. Simonsen, K. A. Eriksen, S. Maslov and K. Sneppen, Physica A 336, 163 (2004).
- (37) M. Cross and H. Greenside, Pattern Formation and Dynamics in Nonequilibrium Systems, (Cambridge, 2009).
- (38) P. N. McGraw and M. Menzinger, Phys. Rev. E 77, 031102 (2008).
- (39) C. A. Lugo and A. J. McKane, Phys. Rev. E 78, 051911 (2008).
- (40) P. de Anna, F. Di Patti, D. Fanelli, A. J. McKane, T. Dauxois, Phys. Rev. E 81, 056110 (2010).