Fluctuation theory in space and time: white noise in reaction-diffusion models of morphogenesis
FFluctuation theory in space and time:white noise in reaction-diffusion models of morphogenesis
Roman Belousov, ∗ Adrian Jacobo, and A. J. Hudspeth
Howard Hughes Medical Institute, Laboratory of Sensory Neuroscience,The Rockefeller University New York, NY 10065, USA (Dated: May 31, 2018; Revision: 1.0)The precision of reaction-diffusion models for mesoscopic physical systems is limited by fluctu-ations. To account for this uncertainty, Van Kampen derived a stochastic Langevin-like reaction-diffusion equation that incorporates spatio-temporal white noise. The resulting solutions, however,have infinite standard deviation.
Ad hoc modifications that address this issue by introducing micro-scopic correlations are inconvenient in many physical contexts of wide interest. We instead estimatethe magnitude of fluctuations by coarse-graining solutions of the Van Kampen equation at a rele-vant mesoscopic scale. The ensuing theory yields fluctuations of finite magnitude. Our approach isdemonstrated for a specific biophysical model—the encoding of positional information. We discussthe properties of the fluctuations and the role played by the macroscopic parameters of the under-lying reaction-diffusion model. The analysis and numerical methods developed here can be appliedin physical problems to predict the magnitude of fluctuations. This general approach can also beextended to other classes of dynamical systems that are described by partial differential equations.
I. INTRODUCTION
In addition to applications in chemistry and other dis-ciplines [1], reaction-diffusion (RD) equations are com-monly accepted as the basis of morphogenetic models inbiology [2–6]. A classical example is the encoding of po-sitional information (PI). During embryological develop-ment, an organism must be partitioned into distinct mor-phological and functional components. The positions ofthese structures may be specified by a chemical agent—a morphogen —whose local concentration varies across theembryo and obeys RD equations. In this context one en-counters perhaps the simplest example of such systems,which has been chosen to demonstrate the theory pre-sented in this paper.As a typical RD system we consider the dynamics ofa single morphogen that diffuses from a localized sourceover a confined spatial domain and undergoes chemicaldegradation (Fig. 1). Once all transients have decayedand the system has reached a steady state, cells or or-ganelles can measure their distance to the source by read-ing out the local concentration of the morphogen. Forthis reason it is said that the morphogen encodes PI.Quantitative characterization of PI noise is an impor-tant problem in biophysics [5, 7–11]. Because many keyprocesses in development occur on micrometer scales, theunderlying chemical reactions and diffusive flows are sub-ject to spontaneous variations. These fluctuations dis-rupt the local concentration of the morphogen and re-duce the amount of information that a RD system con-tains [5, 7, 10]. A relevant question is then: how reliablycan PI be encoded and read out in the presence of noise?Most studies concentrate on the problem of decodingPI. For example, one can estimate the efficiency with ∗ [email protected] Source ( a )0 Lx α ( x ) Distance C on ce n t r a ti on a ( t , x ) FIG. 1. A morphogen produced at the left boundary propa-gates by diffusion into the rest of the one-dimensional systemΛ = [0 , L ]. The morphogen’s steady-state concentration α ( x )owing to a degradation reaction decreases monotonically to-wards the impenetrable right boundary. Each position coor-dinate x ∈ Λ corresponds to a unique value of the positional-information curve α ( x ). The instantaneous concentration ofthe morphogen a ( t, x ), however, is subject to spontaneousfluctuations on a mesoscopic scale. which a cell measures a morphogen’s concentration [7,12]. The concentration can also be measured [7, 8]; theexperimental precision then provides an upper bound forthe uncertainty of PI. Our understanding of the readoutproblem is incomplete, however, for one should also takeinto account how much information a noisy RD systemactually contains.The physical theory of fluctuations opens an avenue tothe problem of encoding PI. At the mesoscopic scale, thedynamics of a reaction-diffusion system can be describedby a stochastic partial differential equation derived fromsimplified microscopic mechanics [13–15]. The noise levelin this model is completely determined by the macro- a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a y scopic parameters of the system, such as the diffusionand reaction constants. The magnitude of fluctuationsin the morphogen’s concentration should then in theorybe calculable. This approach promises clearer results onthe precision of PI than the analysis of empirical data.As shown in Sec. II, the Van Kampen equation leadsto a solution of infinite variance and therefore also ofinfinite standard deviation. Because both of these statis-tics measure the magnitude of fluctuations, one may re-gard this result as futile and seek a more realistic model.The existing alternatives [16–19] either part ways withthe Van Kampen equation or require an additional, adhoc layer of theory. Both approaches, however, rely onnew phenomenological constants such as the amplitude orcorrelation length of microscopic noise. Although theseparameters control and regulate the fluctuation’s mag-nitude, they can be inferred neither from the meso- ormacroscopic dynamics nor from the ensuing theory it-self. Because one can only fit the new parameters toobservations, these models are purely descriptive. Thislack of predictive power is one reason why the theoreticalavenue to the problem of encoding PI has received littleattention.In contradistinction to the theoretical approaches men-tioned above, we estimate the fluctuations in a PI prob-lem by solving the Van Kampen equation without mod-ifications. A plausible level of noise is obtained if theresultant morphogen concentration is integrated in spaceover a subscale of the RD system. This procedure isconsistent with classical fluid dynamics, in which macro-scopic fields are commonly understood as coarse-grainedrepresentations of microscopic systems [20, 21].Multiscale models of computational physics, whichcombine the methods of finite elements and moleculardynamics [22, 23], make the coarse-graining procedureeven more explicit. The macroscopic properties of amolecular-dynamics system are calculated as spatial av-erages by means of a microscopic connection [24, Sec.3 and 6]. The exact procedure amounts to integrationof molecular degrees of freedom over volume, which issuggestive of the coarse-graining subscale. The finite-element method then offers techniques to solve the dy-namical equations of macroscopic fields. Note that, asthe volume of a molecular-dynamics system decreases,the uncertainty of spatial averages diverges, exactly as inthe Van Kampen theory.A coarse-graining subscale arises quite naturally indevelopmental biology: morphogenetic features are notpoint-like, but have a finite mesoscopic extent in space.Moreover, developmental decisions are often delegated towhole cells or to large organelles such as cellular nuclei[7, 11]. In the RD problems of morphogenesis and PI,one should therefore reckon with the total amount of thesubstance and its fluctuations over the scale of the targetbiological structure, rather than with the concentrationfield at isolated points.In the next section we briefly describe the Van Kam-pen equation and the coarse-graining of its solution for a simple RD system. The implications of this theoryand some numerical results are then discussed in Sec. III.Additional mathematical details are provided in the Ap-pendices. In particular, Appendices C and D concern twoclasses of the finite-element method used to simulate nu-merically the dynamics of the coarse-grained stochasticfields. II. THEORY
Because RD problems in general may not yield to an-alytical techniques, we use as a case study our earlier ex-ample of a simple one-dimensional system (Fig. 1). Theassociated theory can be treated by a variety of meth-ods. Purely numerical techniques then can be comparedwith a more accurate analytical approach. This exampleis not entirely abstract, for it provides a model of theactual mechanism of PI encoding in
Drosophila embryos[7, 25].In one dimension the number density of a morphogen a ( t, x ), which depends on time t and position x , obeysthe Van Kampen dynamic equation [13–15]( ∂ t + k − D∂ x ) a ( t, x ) = f ( t, x ) , (1)in which the degradation rate k and the diffusivity D arepositive constants, whereas f ( t, x ) represents microscopicnoise. Appendix A offers a short justification of the VanKampen equation.The left-hand side of Eq. (1) expresses the differencebetween the local change of concentration ∂ t a ( t, x ) andthe classical nonequilibrium forces of mass action andFick’s diffusion. In small systems the residual force f ( t, x ) does not vanish, but varies spontaneously becauseof microscopic events: this is the origin of microscopicnoise. To be consistent with classical fluid dynamics, thesteady-state ensemble averages of f ( t, x ) and a ( t, x ) mustyield (cid:104) f ( t, x ) (cid:105) = 0 , (cid:104) a ( t, x ) (cid:105) = α ( x ) , (2)in which α ( x )—the PI curve—is the time-independentsolution of the macroscopic RD problem [Fig. 1; Eq. (A8)in Appendix A].A convenient model of the morphogen’s source is afixed-value condition imposed at the left end of the in-terval Λ = [0 , L ]. In the macroscopic RD problem, thisconstraint is supplemented quite naturally by a reflectiveright boundary, leading to the expression (A8) for theconcentration curve α ( x ). Nonetheless, in Appendices Aand E we employ a different choice of the right boundarycondition for the stochastic equation (1): a ( t, x ) (cid:12)(cid:12)(cid:12) x =0 = a = α (0) , a ( t, x ) (cid:12)(cid:12)(cid:12) x = L = α ( L ) , (3)with a representing a source of constant strength. Byvirtue of Eq. (2), the fixed-value boundary conditionremains macroscopically consistent with the reflectiveboundary for the PI curve: ∂ x α ( x ) (cid:12)(cid:12)(cid:12) x = L = 0 . (4)The problem of boundary conditions is addressed in Ap-pendix E.At the mesoscale and in the tradition of fluctuationtheory, Van Kampen relates f ( t, x ) to two stochasticterms, owing to the fluctuations of the mass-action lawand the diffusive flow, respectively: f ( t, x ) = (cid:112) kα ( x ) ∂ x ˙ W ( t, x )+ ∂ x [ (cid:112) Dα ( x ) ∂ x ˙ W ( t, x )] . (5)Here ∂ x ˙ W and ∂ x ˙ W are independent, spatially dis-tributed, Gaussian white-noise variates of zero mean andunit strength (Appendix A). The overscript dots indi-cate the time derivatives. These noise sources are delta-correlated in both space and time: (cid:104) ∂ x ˙ W (cid:12)(cid:12)(cid:12) t ,x ∂ x ˙ W (cid:12)(cid:12)(cid:12) t ,x (cid:105) = 0 , (6) (cid:104) ∂ x ˙ W i (cid:12)(cid:12)(cid:12) t ,x ∂ x ˙ W i (cid:12)(cid:12)(cid:12) t ,x (cid:105) = δ ( t − t ) δ ( x − x ) , (7)which hold for i = 1 ,
2, with δ ( · ) being the Dirac deltafunction. Note in the above equations that spatially dis-tributed white noise is singular in time and space: itsvariance diverges as a product of two delta functions,lim t → δ ( t ) and lim x → δ ( x ).To avoid immaterial details, we focus on the steady-state solution of Eq. (1), a ( ∞ , x ). We denote the devia-tion of the morphogen’s concentration from the ensembleaverage value by ∆ a ( t, x ) = a ( t, x ) − α ( x ). Then Eq. (B7)of Appendix B gives us∆ a ( ∞ , x ) = lim t →∞ (cid:90) t dt (cid:48) (cid:90) Λ dx (cid:48) g ( t − t (cid:48) , x | x (cid:48) ) f ( t (cid:48) , x (cid:48) ) . (8)Here g ( t − t (cid:48) , x | x (cid:48) ) is the Green’s function that propagatesdisturbances of the number density in time and space,from an instant t (cid:48) and position x (cid:48) to any other t and x .The steady-state variance of the deviation ∆ a ( ∞ , x ),as can be formally calculated from Eq. (8), diverges [seealso Eq. (B8) in Appendix B]. To understand why thishappens, apply the differential chain rule to the secondterm on the right-hand side of Eq. (5) and substitute itinto Eq. (8); one then finds the following term in theexpression for ∆ a ( ∞ , x ):lim t →∞ (cid:90) t dt (cid:48) (cid:90) Λ dx (cid:48) g ( t − t (cid:48) , x | x (cid:48) ) (cid:112) Dα ( x (cid:48) ) ∂ x (cid:48) ˙ W ( t (cid:48) , x (cid:48) ) ∼ ∂ x W ( t, x ) . (9)Here the time integration removes the temporal singular-ity of ∂ x ˙ W , but the spatial integral is canceled by one ofthe two derivative operators. The above term contains a spatial singularity of the order ∂ x W ( t, x ) [Eq. (7)].Therefore the variance of ∆ a ( ∞ , x ), expressed formallyby Eq. (B8), in effect diverges.If one replaces the positional delta function δ ( · ) inEq. (7) by some bounded correlation kernel C ( · ) [16, Sec.2.1.2], the spatial singularity disappears from Eqs. (8)and (9). The cost of this approach is a significantly morecomplicated theory [26–28]. First, spatial noise correla-tions that regularize the variance of ∆ a ( ∞ , x ) must bemodeled explicitly. Second, a nontrivial kernel C ( · ) in-creases the mathematical difficulty of the problem. Thespatial singularity of Eq. (9) can alternatively be removedby integrating it with respect to the coordinate x , the ap-proach that we pursue here. An integration with respectto position occurs when we coarse-grain the number den-sity a ( t, x ) over a scale ξ of the appropriate spatial di-mension. Then, instead of the morphogen’s concentra-tion at some point x , the quantity of interest becomesthe total number of molecules in the ξ -neighborhoodΞ( x ) = ( x − ξ/ , x + ξ/
2) of that point. If we use theinverse scale ξ as a normalization factor, we can equiva-lently consider a coarse-grained concentration field a ξ ( x ) = (cid:90) x + ξ/ x − ξ/ dx (cid:48) ξ a ( ∞ , x (cid:48) ) . (10)The coarse-grained concentration of the morphogenundergoes fluctuations of finite magnitude. As a sta-tistical measure of this magnitude one can take eitherthe variance of a ξ —the second cumulant κ ( a ξ )—or thestandard deviation std( a ξ ). By using the properties ofthe Green’s and delta functions together with Eqs. (3)–(8), we find κ [ a ξ ( x )] = (cid:90) (cid:90) Ξ( x ) dx dx ξ lim t →∞ (cid:90) t dt (cid:48) (cid:90) Λ dx (cid:48) α ( x (cid:48) ) × (cid:104) kg ( t − t (cid:48) , x | x (cid:48) ) g ( t − t (cid:48) , x | x (cid:48) )+ 2 D∂ x (cid:48) g ( t − t (cid:48) , x | x (cid:48) ) ∂ x (cid:48) g ( t − t (cid:48) , x | x (cid:48) ) (cid:105) . (11)A Fourier series expansion of the above expression, aswell as of the coarse-grained steady-state concentration α ξ ( x ), is derived in Appendix B [Eqs. (B10) and (B11)].The number density α ξ ( x ) differs from α ( x ) by a factorthat is negligible for small scales ξ . Both these fieldsinterchangeably represent a PI curve, for they conveynearly the same value everywhere in Λ.A useful way to quantify the uncertainty of a ξ ( x ) is thecoefficient of variation, std( a ξ ) /α ξ , which relates the levelof fluctuations to the strength of the PI signal. Quite gen-erally, however, both the mean value of the morphogen’sconcentration and its variance are proportional to the pa-rameter a [Eq. (B10) and (B11), Appendix B]. Thereforethe relative uncertainty is inversely proportional to √ a :std[ a ξ ( x )] α ξ ( x ) = σ ξ ( x | (cid:96) ) √ a , (12) (a) (b) FIG. 2. Uncertainty of the coarse-grained concentration α ξ ( x ) as a function of position for the experimental data of Refs. [7, 25].The values of the standard deviation in A and the coefficient of variation in B are calculated by three methods: i) analyticalsolution (Appendix B); ii) spectral finite-element simulations (Appendix C); iii) collocation finite-element method (Appendix D).The uncertainty of the PI curve does not exceed 0 . in which σ ξ ( x | (cid:96) ) depends only on the coordinate x andthe parameters ξ and (cid:96) = (cid:112) D/kL − .Fluctuations of physical quantities usually decay as theinverse square root of the number of molecules involved[29]. This dependence is explicitly controlled by the pa-rameter a in Eq. (12). The source strength a in theabove expression should therefore be measured in one di-mension as the number of molecules per unit length. Ifmolar or mass-density units are used instead, Eq. (12)does not render the coefficient of variation correctly.Equation (12) defines σ ξ ( x | (cid:96) ), which we call a varia-tion profile . Given the values of ξ and (cid:96) , this relationcan be evaluated numerically as a function of position byuse of Eqs. (B10) and (B11). For a source of any givenstrength a , the coefficient of variation—a rescaled varia-tion profile—can then be calculated easily from Eq. (12).Finally, the constant (cid:96) > a ξ ( x ) [Appendix B, Eq. (B13)]. The fluctuations ofthe morphogen’s concentration at two points separatedby distances greater than (cid:96)L = (cid:112) D/k are nearly inde-pendent, whereas the decay of the time correlations iscontrolled by ( k(cid:96) ) − = L /D (Appendix B). The con-stants k and L determine the scale of the system; theycan serve as units of time and length. III. NUMERICAL RESULTS
As an application of the Van Kampen theory, we esti-mate the level of fluctuations for the concentration of themorphogen bicoid in a
Drosophila embryo [7, 25]. The re-sults of our calculations are reported in a system of unitsreduced by the length constant L and the time constant k − . The values of the physical parameters are adoptedfrom experimental data [7, 25]: L = 0 . λ = 0 . a = 55 nM (1 nM corresponds to 0.6 molecules/ µ m ).Converted to reduced units, the source strength and thecorrelation length are respectively a = 4 . × and (cid:96) = 0 .
2. The concentration of bicoid is presumably readout by densely distributed cellular nuclei, whose spatialseparation sets a plausible coarse-grain scale of ξ = 0 . a ξ ) on the position x , calculatedfor the coarse-grained PI curve. This is a convex curvethat is defined over the interval [ ξ/ , L − ξ/
2] and de-creases monotonically from its maximum at x = ξ/ λ = (cid:96)L , the variance of thecoarse-grained morphogen concentration depends almostlinearly on x and flattens when λ → ∞ and k →
0. Inthe latter case, which represents pure diffusion withoutdegradation, the fluctuations of α ( x ) are maximal for anygiven values of a > ξ >
0. On the other hand,when the rate constant k becomes infinitely large, theproblem degenerates and fluctuations vanish.Of the two numerical integration schemes discussedin this paper, the collocation method (Appendix D) isless accurate than the spectral finite-element algorithm(Appendix C). The latter approach compares favorablywith the analytical solution [Eq. (B10) in Appendix B].Both integration algorithms nonetheless reproduce cor-rectly the average concentration and the overall shape ofthe variation profile.The relative uncertainty of the PI curve in our nu-merical example does not exceed 0 . α ξ decreases with x faster than itsstandard deviation, the coefficient of variation increasestowards the right boundary. The relative uncertaintynevertheless remains within the order of 0 . k(cid:96) ) − . This result validates the Van Kam-pen theory for conditions approaching the macroscopiclimit. However, in developmental processes on a scale oftens of micrometers [30]—the dimension relevant to thespecification of intracellular structures—fluctuations canchallenge the efficiency of morphogen receptors. Addi-tional mechanisms, such as biochemical feedback loops[31], might then be required to reduce noise in the sys-tem. IV. CONCLUSION
The Van Kampen theory provides a promising meansof estimating the fluctuation level in RD problems andmore generally in systems of mesoscopic physical fields.The approach is conceptually simple and has a relativelysmall computational cost. Although in this paper weconsider only the steady-state solution of a RD problem,transients can be taken into account as well [Eq. (B6) inAppendix B]. Moreover, the Van Kampen theory can beintegrated readily into multiscale computational models.To simulate the Van Kampen equation, we formulatedand tested two numerical techniques. The results of spec-tral finite-element simulations (Appendix C) were quiteaccurate and superior to those of the collocation method(Appendix D).As a case study we chose a relatively large, 500 µ m-long system for its simple geometry and the availability ofexperimental data. Because the length scale approachesmacroscopic conditions, fluctuations of the PI curve inour simple example are very small. In many other in-stances of the PI problem, however, the system’s size canbe 10 µ m or even less. At such scales, the fluctuationsof the PI curve can impose operational time and spaceconstraints on the detectors of morphogen concentration.Estimation of the noise level might provide insight intothe mechanisms of encoding and readout of PI. For exam-ple, the concentration’s uncertainty might help in identi-fying a morphogen among the candidate substances thatoccur in a system.In a study focused on a specific RD problem, there aremore details that could be included in a Van Kampenequation: fluctuations of the source strength, boundaryeffects, and the dimensionality of the system. Incorpo-ration of these factors should improve the accuracy of atheoretical model (Appendices B and E). ACKNOWLEDGMENTS
The authors thank Dr. A. Erzberger, Dr. A. Milewski,and Dr. F. Berger for stimulating discussions and valu-able comments on our results.A. J. is a Fellow of the F. M. Kirby Foundation. R. B.is a Research Associate and A. J. H. an Investigator of Howard Hughes Medical Institute.
Appendix A: Van Kampen reaction-diffusionequation
The Van Kampen RD equation (1) extends theLangevin model of fluctuations for simple time-dependent physical quantities to spatially distributedfields [13–15]. Consider first the classical RD dynamicsfor the number density a ( t, x ) of some morphogen overthe linear domain x ∈ Λ: ∂ t a ( t, x ) = − ka ( t, x ) + D∂ x a ( t, x ) , (A1)in which the degradation rate k and the diffusivity D arepositive constants. The first term on the right-hand sideof Eq. (A1) states the mass-action law for the chemicaldegradation of the morphogen. The second term repre-sents the divergence of the Fick’s diffusion flow J ( t, x ) = − D∂ x a ⇒ − ∂ x J ( t, x ) = D∂ x a, (A2)which describes the balance of incoming and outgoingcurrents of matter J ( t, x ).Both macroscopic laws—mass action and Fick’sdiffusion—emerge as statistical averages of the micro-scopic dynamics [13, 14] over a steady-state ensemble.Ergodicity of a system is commonly assumed as well. Atmesoscopic scales, however, we must allow fluctuationsby replacing the following terms in Eqs. (A1) and (A2): ka ( t, x ) → ka ( t, x ) − χ ( t, x ) ,J ( t, x ) → − D∂ x a ( t, x ) − j ( t, x ) , (A3)in which χ ( t, x ) and j ( t, x ) are the deviations from theclassical macroscopic laws of reaction and diffusion, re-spectively. Due to the spontaneous variations given byEqs. (A3), the local change of concentration ∂ t a ( t, x ) doesnot exactly match the fluctuating force on the right-handside of Eq. (A1). The residual is( ∂ t + k − D∂ x ) a ( t, x ) = χ ( t, x ) + ∂ x j ( t, x ) . (A4)The spontaneous behavior of the fluctuating force on theright-hand side of this equation appears practically ran-dom and is therefore modeled as a stochastic, spatiallydistributed process.In the Langevin approach, the macroscopic propertiesof a steady-state dynamics correspond to the ensembleaverages, here denoted by angle brackets, of mesoscopicvariables. In particular, Eq. (A1) requires (cid:104) χ ( t, x ) (cid:105) = 0and (cid:104) ∂ x j ( t, x ) (cid:105) = 0. Hence the ensemble average ofEq. (A4) yields k (cid:104) a ( t, x ) (cid:105) − D∂ x (cid:104) a ( t, x ) (cid:105) = 0 (A5)for (cid:104) ∂ t a ( t, x ) (cid:105) = 0 by the definition of a steady state.For a nontrivial solution (cid:104) a ( t, x ) (cid:105) (cid:54) = 0 to exist, wemodel a source of the chemical agent by a nonhomoge-neous Dirichlet boundary condition, which fixes the valueof (cid:104) a ( t, x ) (cid:105) at x = 0: (cid:104) a ( t, (cid:105) = a . (A6)A natural choice of the other boundary at x = L is thehomogeneous Neumann condition, which reflects the dif-fusive flow (cid:104) J ( t, x ) (cid:105) [Eq. (A3)]: (cid:104) ∂ x a ( t, L ) (cid:105) = 0 . (A7)Subject to the above constraints, Eq. (A5) is easy tosolve [32, Chapter 2]. We thus obtain the macroscopictime-independent PI curve plotted in Fig. 1, α ( x ) = (cid:104) a ( t, x ) (cid:105) = a cosh[( L − x ) /λ ]cosh( L/λ ) , (A8)in which λ = (cid:112) D/k .For Langevin dynamics (A7) we reduce the homoge-neous Neumann condition at the right end to a consistentnonhomogeneous Dirichlet boundary: (cid:104) a ( t, L ) (cid:105) = α ( L ) . (A9)As discussed in Appendix E, both Dirichlet condi-tions (A6) and (A9) neglect fluctuation effects at theboundaries of the RD system. These effects can be in-cluded by imposing the reflective Neumann conditions onboth ends of Λ. These details, which are not strictly nec-essary for a simple demonstration of the Van Kampentheory, are spared for Appendix E.The Langevin model is complete once the statisticalproperties for the right-hand side of Eq. (A4) have beenspecified. Van Kampen derives them by reducing themicroscopic RD dynamics to a random walk, a tradi-tional argument of statistical mechanics [33]. A continu-ous limit of this simplified model gives (cid:104) χ ( t, x ) χ ( t (cid:48) , x (cid:48) ) (cid:105) = kα ( x ) δ ( t − t (cid:48) ) δ ( x − x (cid:48) ) , (A10) (cid:104) j ( t, x ) j ( t (cid:48) , x (cid:48) ) (cid:105) = 2 Dα ( x ) δ ( t − t (cid:48) ) δ ( x − x (cid:48) ) , (A11)which hold for any instants of time t, t (cid:48) , and positions x, x (cid:48) [13, 14]. The theory behind the above equations re-lies on the following assumptions: χ and j are indepen-dent ( (cid:104) χ ( t, x ) j ( t (cid:48) , x (cid:48) ) (cid:105) ≡ dx contains a large number of molecules α ( x ) dx ; and all cor-relations at distances of order dx are negligible. Then in-finitesimal processes χ ( t, x ) and j ( t, x ) are approximatelyGaussian by virtue of the central limit theorem [34, Sec.2.5].The Van Kampen model leads directly to the conceptof a spatially distributed Gaussian white noise ∂ x ˙ W ( t, x )with a zero mean and a constant strength β . The definingproperty of ∂ x ˙ W is that its integral over a time interval t and a line segment Ξ( x ) = ( x − ξ/ , x + ξ/ W ( t, x | ξ ) = (cid:90) t dt (cid:48) (cid:90) Ξ( x ) dx (cid:48) ∂ x (cid:48) ˙ W ( t (cid:48) , x (cid:48) ) , (A12) is a Gaussian random process of zero mean ( (cid:104) W (cid:105) = 0)and variance (cid:104) W ( t, x | ξ ) (cid:105) = βξt. (A13)All properties of the stochastic processes χ ( t, x ) and j ( t, x ) are then encompassed by χ ( t, x ) = (cid:112) kα ( x ) ∂ x ˙ W ( t, x ) , (A14) j ( t, x ) = (cid:112) Dα ( x ) ∂ x ˙ W ( t, x ) , (A15)in which ∂ x ˙ W and ∂ x ˙ W are two independent, spa-tially distributed, Gaussian white-noise terms of unitstrength β = 1 [Eqs. (6) and (7)]. Note that j ( t, x ) isa vector quantity, which in one dimension has a singlecomponent ∂ x ˙ W ( t, x ). Appendix B: Green’s function method
Supplemented with an initial value a (0 , x ) and theboundary conditions (3), Eqs. (A4)–(A15) lead toEq. (1): ( ∂ t + k − D∂ x ) a ( t, x ) = f ( t, x ) . This stochastic partial differential equation is linear, asis its left-hand-side operator L = ( ∂ t + k − D∂ x ), andinhomogeneous, in that f ( t, x ) enters the expression ad-ditively.A general solution of Eq. (1) is most conveniently ex-pressed with the help of the Green’s function [35, Chap-ter 10] g ( t − t (cid:48) , x | x (cid:48) ), which we find from the followingequations: L g ( t − t (cid:48) , x | x (cid:48) ) = δ ( t − t (cid:48) ) δ ( x − x (cid:48) ) , (B1) g ( t,
0) = g ( t, L ) = 0 . (B2)In a finite domain Λ the Green’s function can be ex-panded as a series [36]. For the problem at hand weuse a discrete expansion ( n = 1 , ... ) in the orthonormalFourier basis φ n ( x ) = (cid:112) /L sin( nπx/L ) , (B3)which is complete under the boundary conditions (B2).One then finds g ( t − t (cid:48) , x | x (cid:48) ) = ∞ (cid:88) n =1 g n ( t − t (cid:48) ) φ n ( x (cid:48) ) φ n ( x ) , (B4) g n ( t ) = H ( t ) exp (cid:26) − t (cid:20) k + D n π L (cid:21)(cid:27) , (B5)in which H ( · ) stands for the Heaviside step function.With the boundary conditions (3) and (B2), the gen-eral solution of (1) takes the form a ( t, x ) = α ( x ) + (cid:90) Λ dx (cid:48) g ( t, x | x (cid:48) )[ a (0 , x (cid:48) ) − α ( x (cid:48) )]+ (cid:90) t dt (cid:48) (cid:90) Λ dx (cid:48) g ( t − t (cid:48) , x | x (cid:48) ) f ( t (cid:48) , x (cid:48) ) , (B6)in which the second term on the right-hand side van-ishes as lim t →∞ g ( t, x | x (cid:48) ) = 0. If one is concerned merelywith the steady-state behavior of Eq. (1), the transientsolutions can be eliminated from Eq. (B6) by taking thelimit of infinite t : a ( ∞ , x ) = α ( x )+ lim t →∞ (cid:90) t dt (cid:48) (cid:90) Λ dx (cid:48) g ( t − t (cid:48) , x | x (cid:48) ) f ( t (cid:48) , x (cid:48) ) . (B7)Consider the statistical properties of the steady-statesolution a ( ∞ , x ). Because f ( t, x ) given by Eqs. (5) is alinear superposition of zero-mean, Gaussian white-noiseterms, the ensemble average of Eq. (B7) is consistent withthe macroscopic dynamics (A1): (cid:104) a ( ∞ , x ) (cid:105) = α ( x ) . The second cumulant κ [ a ( ∞ , x )] of a ( ∞ , x ) can be ob-tained from Eqs. (3)–(8), and (B2): κ [ a ( ∞ , x )] = κ [∆ a ( ∞ , x )] = (cid:104) a ( ∞ , x ) − α ( x ) (cid:105) = k lim t →∞ (cid:90) t dt (cid:48) (cid:90) Λ dx (cid:48) α ( x (cid:48) )[ g ( t − t (cid:48) , x | x (cid:48) )] + 2 D lim t →∞ (cid:90) t dt (cid:48) (cid:90) Λ dx (cid:48) α ( x (cid:48) )[ ∂ x (cid:48) g ( t − t (cid:48) , x | x (cid:48) )] . (B8)Higher-order cumulants of the steady-state solution van-ish due to the Gaussian nature of f ( t, x ) and hence of α ( ∞ , x ) as well.As explained in Sec. II, the formal expression (B8) di-verges. Therefore, to estimate the magnitude of fluc-tuations in the morphogen’s concentration, we calcu-late the variance of the coarse-grained number density a ξ ( ∞ , x ) from Eq. (11). This computation can be carriedout through a series expansion of the Green’s function, g ( t − t (cid:48) , x | x (cid:48) ), truncated at N th term. Let us introducethe following formulas:Φ n ( x ) = (cid:90) Ξ( x ) dx (cid:48) ξ φ n ( x (cid:48) ) = 2 Lnπξ sin (cid:18) nπξ L (cid:19) φ n ( x );Ω mn = 4 π (cid:96) mn tanh[ (cid:96) − ][1 + π (cid:96) ( m − n ) ][1 + π (cid:96) ( m + n ) ] , (B9)in which (cid:96) = (cid:112) D/kL − . Then, by substitutingEqs. (A8), (B4), and (B5) into (11) and completing theintegrals, we obtain κ [ a ξ ( x )] = a (cid:88) mn Ω mn Φ m ( x )Φ n ( x ) , (B10)in which the summation runs over all positive integers m and n up to N ( m, n = 1 ..N ). Note that the mean valueof the coarse-grained field a ξ ( x ) is (cid:104) a ξ ( x ) (cid:105) = 2 λξ sinh (cid:18) ξ λ (cid:19) α ( x ) → ξ → α ( x ) . (B11) We can similarly obtain the autocorrelation function κ [ a ξ (0 , x ) , a ξ ( t, x )] for the time-dependent concentra-tion a ξ ( t, x ) = (cid:90) Ξ( x ) dx (cid:48) a ( t, x (cid:48) ) = α ξ ( x ) + ∆ a ξ ( t, x ) . (B12)In linear systems the decay of the temporal and spatialautocorrelations is encompassed by the Green’s function[37, Sec. 8.6]: κ [ a ξ (0 , x ) , a ξ ( t, x )] = (cid:104) ∆ a ξ (0 , x )∆ a ξ ( t, x ) (cid:105) = (cid:90) Λ dx (cid:48) g ( t, x | x (cid:48) ) (cid:104) ∆ a ξ (0 , x )∆ a (0 , x (cid:48) ) (cid:105) = a (cid:88) mn Ω mn Φ m ( x )Φ n ( x ) exp[ − kt (1 + π (cid:96) m )] , (B13)in which only the transient term of Eq. (B6) makes a non-zero contribution [38]. From Eqs. (B4), (B5), (B9), and(B13) it follows that, in reduced units (Sec. III), the timeand space correlations are controlled respectively by theparameters ( k(cid:96) ) = D/L and λ = (cid:96)L = (cid:112) D/k throughthe diffusion constant D .In computations the series expansion (B10) and (B13)should be truncated at an order N ≥ (cid:100) L/ξ (cid:101) . This opti-mal value is suggested by the following argument. Sup-pose that 2
L/ξ is an integer. The Fourier mode φ n + N ( · )is then an alias of φ n ( · ), because φ n + N ( x ) = φ n ( x ) when-ever x is an integer multiple of ξ . In Eq. (B10) we passedfrom the basis set φ n ( · ) to the coarse-grained functionsΦ n ( · ) by integrating over a spatial scale ξ [Eq. (B9)].This procedure allows us to disregard aliasing modeswith n > N . Indeed, spatial features smaller than thescale ξ should be smoothed by the coarse-graining in-tegration. The regions near the ends of the domain Λ( x ≈ ξ/ , L − ξ/
2) are exceptions that may require moreterms to reduce ringing artifacts.
Appendix C: Spectral method
Modal analysis similar to that of Appendix B leadsto a simple method of spectral finite elements [39–41]for Eq. (1). Subject to the boundary conditions (3),the number density a ( t, x ) has a series representation interms of the basis functions given by Eq. (B3): a ( t, x ) = α ( x ) + ∞ (cid:88) m =1 a m ( t ) φ m ( x ) , (C1)in which the time-dependent coefficients a m ( t ) must van-ish on average to satisfy Eq. (2): (cid:104) a m ( t ) (cid:105) = 0.Spatially distributed white noise ∂ x ˙ W i ( i = 1 ,
2) like-wise has a representation ∂ x ˙ W i ( t, x ) = ∞ (cid:88) n =1 ˙ w in ( t ) φ n ( x ) , (C2)in which each time-dependent coefficient ˙ w in ( t ) is a sim-ple, independent Gaussian white noise. Equations (6),(7), (A12), and (A13) readily follow from Eq. (C2). Inhigher dimensions there are additional vector compo-nents like ∂ x ˙ W (Appendix B) that are independent inan orthogonal reference frame and therefore can be ex-panded separately in series (C2). In nonorthogonal co-ordinate systems one must also account for correlationsdue to the overlap of basis vectors.Equations (C1) and (C2), truncated at some m and n ,provide the finite-element representations of a ( t, x ) andthe white-noise terms. Their substitution into (1) and afew simple manipulations eventually lead to a system ofordinary differential equations for the coefficients a m ( t ):˙ a m ( t ) = − k (1 + π (cid:96) m ) a m ( t ) + f m ( t ) , (C3)in which (cid:96) = (cid:112) D/kL − , and f m ( t ) = (cid:88) n [ F mn ˙ w n ( t ) + F mn ˙ w n ( t )] , (C4) F mn = (cid:90) Λ dx (cid:112) kα ( x ) φ m ( x ) φ n ( x ) , (C5) F mn = − (cid:90) Λ dx (cid:112) Dα ( x ) ∂ x φ m ( x ) φ n ( x ) . (C6)From Eq. (C3) we see that stochastic forces f m ( t ) ran-domly perturb the modal coefficients a m ( t ). When ex-plicit analytical expressions are not available for the spa-tial integrals in Eqs. (C5) and (C6), a discrete Fouriertransform can be used instead as an approximation. Thisapproach should then be termed a pseudospectral finite-element method.For a general RD system, the derivation of equationsanalogous to (C3) is quite straightforward. Because theproblem studied in this paper is relatively simple, wecan obtain each coefficient a m ( t ) explicitly (Appendix B).Because numerical methods are more widely applicablethan analytical ones, however, we develop below a pseu-dospectral finite-element scheme to solve the system ofequations (C3).The system of equations (C3) can be numerically in-tegrated in time by various methods, such as the Crank-Nicolson algorithm discussed in the next section of theAppendix. We can also use a second-order stochasticoperator-splitting technique [42, Appendix C]: a m ( t + ∆ t ) = exp[ − (1 + π (cid:96) m ) k ∆ t ] a m ( t )+ exp (cid:20) − (1 + π (cid:96) m ) k ∆ t (cid:21) F m (∆ t ); F m (∆ t ) = (cid:90) ∆ t dt f m ( t )= √ ∆ t (cid:88) n [ F mn w n + F mn w n ] , (C7)in which w n and w n are independent Gaussian randomvariables of zero mean and unit variance, whereas ∆ t isan integration time step. By coarse-graining Eq. (C1), one easily finds a finite-element representation of ∆ a ξ ( x ) = a ξ ( x ) − α ξ ( x ) in thenotation of Eq. (B9):∆ a ξ ( x ) = (cid:88) m a m Φ m ( x ) , (C8)in which the coefficients a m are correlated Gaussian vari-ables. Their covariance matrix K mn = (cid:104) a m a n (cid:105) can be sampled in a numerical simulation of Eq. (C7).The variance of ∆ a ξ ( x ) is then equal to κ [∆ a ξ ( x )] = (cid:88) mn K mn Φ m ( x )Φ n ( x ) . (C9)By comparing the above equation with (B10), we see that K mn = a Ω mn .The above results show that the modal coefficients a m ( t ) are correlated Gaussian random variables of fi-nite mean and variance. However, there are so manyof them in the representation of a ( ∞ , x ) that its variancegiven by Eq. (B8) diverges unless the coarse-grained ba-sis functions Φ m ( x ) are used as in (B10). Thanks to thenonsingular nature of the modal coefficients, a computersimulation of Eq. (C7) is feasible.The series expansions (C1) and (C2) need not have thesame number of modes. The argument of Appendix Bfor Eq. (B10) sets the optimum for m = 1 , , ..N m , inwhich N m = (cid:100) L/ξ (cid:101) . The accuracy of simulations can beimproved indefinitely by taking progressively more termsin Eq. (C2), N n ≥ N m . In Sec. III we report our resultsfor N n = 4 N m and ∆ t = 10 − . Appendix D: Collocation method
The Fourier components a m ( t ) of the morphogen’s con-centration are Gaussian random variables of finite vari-ance (Appendix C). This vector representation of thefield a ( t, x ) can be projected onto another basis set.Then, in principle, it should be possible to simulate nu-merically the dynamics of the new components.In this section we use a piecewise-linear interpola-tion [39, Chapter 1], as a basis set for the finite-volumemethod, a widely used collocation finite-element scheme[43, 44, Chapter 4]. Consider a uniform grid x i = i ∆ x, i = 0 , , . . . M + 1 on the domain Λ (Fig. 3). Wecenter control elements of size ∆ x at the nodes x i , i =1 , , . . . M . The morphogen’s concentration is then inter-polated by a ( t, x ) ≈ M +1 (cid:88) i =0 A i ( t ) η i ( x ) , (D1) η i ( x ) = (cid:40) ∆ x − x i + x ∆ x , if x i − ≤ x ≤ x i ∆ x + x i − x ∆ x , if x i ≤ x ≤ x i +1 , (D2) i =0 i = M +1 i=2 Control element A ( x i ) a xA ( x ) α ( L ) FIG. 3. A piecewise-linear interpolant A ( x ) on a grid x i , i = 0 ..M + 1. We use the boundary conditions (3) tofix the values A ( x ) = a and A ( x M +1 ) = α ( L ) at the nodes i = 0 and i = M + 1 ( ◦ ), respectively. The control elementsare centered on the nodes i = 1 ..M ( • ). This collocationscheme reserves small intervals near the boundaries of the do-main for a convenient implementation of the coarse-grainingprocedure. which coincides with a ( t, x ) at the centers of the controlelements A i ( t ) = a ( t, x i ). To comply with the bound-ary conditions (3), we fix A = a and A M +1 = α ( L ).Thus the above interpolation is completely specified by M time-dependent components A i ( t ) , i = 1 ..M .A standard procedure of the finite-volume methodwould be to integrate Eq. (1) over each i th control el-ement X i = [ x i − ∆ x/ , x i + ∆ x/ ξ ≥ ∆ x inorder to remove the spatial singularity of the stochasticnoise (Sec. II). It is convenient to partition the domainΛ so that ξ is an integer multiple of ∆ x : ξ = P ∆ x . Weconstruct an integral operator I P i = P − P − (cid:88) j =0 I i + j = P − P − (cid:88) j =0 (cid:90) (cid:90) X i + j dx ∆ x . (D3)In the above expression it suffices to consider only P = 1with a single term I i that can be used to evaluateEq. (D3) for an arbitrary P . Applied to a ( t, x ), the op-erator I i gives a spatial integral of the coarse-grainedconcentration: I i a ( t, x ) = (cid:90) X i dx (cid:90) x +∆ x/ x − ∆ x/ dx ∆ x a ( t, x )= (cid:90) X i dx a ξ ( t, x ) (cid:12)(cid:12)(cid:12) ξ =∆ x . (D4)By applying each of the M operators I i , i = 1 , . . . M to both sides of Eq. (1), we get( ∂ t + k ) I i a ( t, x ) − D ∆ x [ a ( t, x i − ∆ x ) − a ( t, x i ) + a ( t, x i + ∆ x )]= I i f ( t, x ) . (D5) Then substituting Eq. (D1) for a ( t, x ) yields∆ x ( ∂ t + k ) (cid:20) A i ( t )3 + A i − ( t ) + A i +1 ( t )6 (cid:21) − D ∆ x [ A i − ( t ) − A i ( t ) + A i +1 ( t )]= I i f ( t, x ) , (D6)which forms a system of M ordinary differential equa-tions to be solved for A i ( t ) , i = 1 , . . . M .The statistical properties of the coarse-grained field a ξ ( t, x ) can be estimated from the values A i ( t ) [Eq. (D1)]sampled in a computer simulation. For example, when ξ = ∆ x one obtains a ∆ x ( t, x i ) ≈ A i ( t )4 + A i − ( t ) + A i +1 ( t )8 ,a ∆ x ( t, x i + ∆ x/ ≈ A i ( t )2 + A i +1 ( t )2 . (D7)To integrate Eq. (D6) in time we use the Crank-Nicolsonscheme based on the trapezoidal rule (cid:90) t +∆ tt dtA i ( t ) ≈ [ A i ( t ) + A i ( t + ∆ t )] ∆ t . (D8)The simulation algorithm can be concisely written in thevector-matrix notation T A ( t + ∆ t ) = EA ( t ) + R (∆ t ); (D9) T ii = (2 + k ∆ t ) ∆ x D ∆ t ∆ x , (D10) T i ( i ± = (2 + k ∆ t ) ∆ x − D ∆ t x , (D11) E ii = (2 − k ∆ t ) ∆ x − D ∆ t ∆ x , (D12) E i ( i ± = (2 − k ∆ t ) ∆ x
12 + D ∆ t x , (D13)with all other elements of M -by- M matrices T and E being zero. The column vector R (∆ t ) of size M is givenbelow.On the right-hand side of Eq. (D6) we obtain two con-tributions, which correspond to the stochastic terms ofthe force f ( t, x ) [Eq. (5)]: (cid:90) ∆ t dt I i (cid:112) kα ( x ) ∂ x ˙ W ( t, x ) = u i , (cid:90) ∆ t dt I i ∂ x [ (cid:112) Dα ( x ) ∂ x ˙ W ( t, x )] = v i − v i +1 , (D14)in which u i and v i are Gaussian random variables of zeromean. The covariances of u i and v i can be calculatedfrom Eq. (6) and (7). In particular, all v i are independentof each other, as well as from u i which correlate in pairs: (cid:104) u i u j (cid:105) (cid:54) = 0 if | i − j | = 1. We omit lengthy calculations0that eventually lead to κ ( u i ) = 4 λ k ∆ tα ( x i )∆ x (cid:20) sinh (cid:18) ∆ xλ (cid:19) − ∆ xλ (cid:21) , (D15) (cid:104) u i u i +1 (cid:105) = 4 λ k ∆ ta ∆ x cosh( L/λ ) cosh (cid:18) L − x i − x i +1 λ (cid:19) × (cid:20) ∆ x λ cosh (cid:18) ∆ x λ (cid:19) − sinh (cid:18) ∆ x λ (cid:19)(cid:21) , (D16) κ ( v i ) = 2 λD ∆ ta ∆ x cosh( L/λ ) × (cid:20) sinh (cid:18) L − x i − λ (cid:19) − sinh (cid:18) L − x i λ (cid:19)(cid:21) , (D17)with λ = (cid:112) D/k .An alternative approach can be used for prob-lems in which the derivation of formulas analogous toEqs. (D15)–(D16) is too tedious. If we substitutethe spectral representation of white noise (C2) intoEq. (D14), we get u i = √ ∆ t (cid:88) j w j u ij , v i = √ ∆ t (cid:88) j w j v ij ,u ij = I i [ (cid:112) kα ( x ) φ j ( x )] ,v ij = (cid:90) X i dx ∆ x (cid:112) Dα ( x + ∆ x/ φ j ( x + ∆ x/ , (D18)in which w j and w j are independent Gaussian randomvariables with zero mean and unit variance. The ele-ments u ij and v ij can be evaluated by the pseudospectralmethod (Appendix C). In essence, the above equationsare projections of white noise in a spectral representationonto the linear-interpolation basis functions, a procedurementioned in the beginning of this section.Finally, we can express the column vector R as R = b + u + Cv , (D19) b = a ∆ t (cid:18) D ∆ x − k ∆ x (cid:19) , (D20) b M = α ( L )∆ t (cid:18) D ∆ x − k ∆ x (cid:19) − v M +1 (D21) C ii = 1 , C i ( i +1) = − , (D22)in which the vector b incorporates the boundary terms;the components b ij and the matrix elements C ij thatequal zero are not indicated.The results reported in Sec. III are obtained by usingEqs. (D15)–(D17) with ∆ t = 10 − and ∆ x = 2 × − . Appendix E: Fluctuation effects at the source andboundaries
The Dirichlet conditions (3) fix the value of a ( t, x )at the ends of the domain Λ. The resulting solution of Eq. (1) therefore neglects fluctuation effects at theboundaries and, in particular, at the source of the mor-phogen. Indeed, the Green’s function we obtained in Ap-pendix B is insusceptible to any forces at the ends of thedomain Λ, where it vanishes due to Eq. (B2). Note thatthe left Dirichlet boundary in Eq. (3) is an implicit sourceof the morphogen.Assuming a closed RD system, in which matter doesnot leak through the ends of the domain Λ, we replacethe reflective Neumann boundary conditions (3) by ∂ x a ( t, x ) (cid:12)(cid:12)(cid:12) x =0 = 0 , ∂ x a ( t, x ) (cid:12)(cid:12)(cid:12) x = L = 0 , (E1)whereby the macroscopic diffusion flux through thepoints x = 0 , L vanishes [Eq. (A2)]. Nonetheless we shalltake special care that the fluctuations of the matter flow j ( t, x ) in Eq. (A3) do not violate the closed-system con-straint.Together with Eq. (E1), we must model explicitly thesource of the morphogen s ( x ). Given that this substanceis generated from a densely concentrated substrate at aneffective rate k + , one can pose s ( x ) = k + δ ( x ) + (cid:112) k + ˙ w + ( t ) δ ( x ) , (E2)in which the second term, with simple white noise co-efficient ˙ w + ( t ), introduces fluctuations of the sourcestrength.The new boundary conditions (E1) require a change ofthe basis set φ n ( · ) → ψ n ( · ) in the series expansion for theGreen’s function [Eq. (B4)], for the concentration a ( t, x )[Eq. (C1)], and finally for white noise ∂ x ˙ W , but not for ∂ x ˙ W . The Neumann basis set has an additional modefor n = 0: ψ ( x ) = L − / , ψ n ( x ) = (cid:112) /L cos( nπx/L ) . (E3)The Neumann boundary conditions also alter the expres-sion for the PI curve α ( x ) → ν ( x ): ν ( x ) = k + k ∞ (cid:88) n =0 ψ n (0) ψ n ( x )1 + π λ n /L . (E4)White noise in the expression for the fluctuating diffu-sive flow, ∂ x ˙ W ( t, x ), should be treated separately. Be-cause this term models variations of matter flux througha point x , its spatial derivative at x = 0 , L is not welldefined. With the no-leak conditions, matter is not al-lowed to flow through the boundaries of the domain Λ.To satisfy this requirement we must use the expansionEq. (C2) for ∂ x ˙ W ( t, x ).If we were to enforce the expansion of ∂ x W ( t, x ) inthe basis set (E3) instead of (B3), through integrationby parts in Eq. (8) we would obtain a boundary term ofthe formlim t →∞ (cid:90) t dt (cid:48) g ( t − t (cid:48) , x | x (cid:48) ) (cid:112) Dα ( x (cid:48) ) ∂ x (cid:48) W ( t (cid:48) , x (cid:48) ) (cid:12)(cid:12)(cid:12) Lx (cid:48) =0 . (E5)1As for Eq. (9), this expression has a spatial singular-ity because the Green’s function does not vanish at theNeumann boundaries. We cannot remedy this problemby coarse graining, for we cannot integrate over a ξ -neighborhood of the end points x = 0 , L .If the morphogen is allowed to leak through the bound-aries of the system, new point sources of fluctations, sim-ilar to the second term in Eq. (E2), may emerge. TheVan Kampen theory should then be extended to includethis contribution.In summary, we introduced above three point sourcesof fluctuations: noise in the morphogen’s synthesis at x = 0 and in its degradation at x = 0 and x = L . Inour simplified model of Sec. II, however, there is a con- tinuum of noise over the interval x ∈ (0 , L ). Providedthe average total number of molecules in the system islarge, we do not expect that the addition of three isolatedpoints, as suggested in this section, would alter the levelof fluctuations.Systematically applying the above changes to the the-ory presented earlier, one can derive results that incor-porate fluctuations at the source of the chemical agentand at the boundaries of the RD system. From theanalysis of this section it also follows that the Dirichlet-Neumann conditions (A6) and (A7) would not providea full account of these phenomena. A more reasonablechoice would be either to increase the level of details byEqs. (E1) and (E2) or to neglect the boundary effectsaltogether by using Eq. (3). [1] G. Nicolis and A. D. Wit, Scholarpedia , 1475 (2007),revision , 37 (1952).[3] L. Wolpert, Journal of Theoretical Biology , 1 (1969);in Essays on Developmental Biology, Part B , CurrentTopics in Developmental Biology, Vol. 117, edited byP. M. Wassarman (Academic Press, 2016) pp. 597 – 608.[4] J. B. A. Green and J. Sharpe, Development , 1203(2015).[5] G. Tkaˇcik, J. O. Dubuis, M. D. Petkova, and T. Gregor,Genetics , 39 (2015).[6] J. Halatek and E. Frey, Nature Physics , 507 (2018).[7] T. Gregor, D. W. Tank, E. F. Wieschaus, and W. Bialek,Cell , 153 (2007).[8] T. Gregor, E. F. Wieschaus, A. P. McGregor, W. Bialek,and D. W. Tank, Cell , 141 (2007).[9] T. Bollenbach, P. Pantazis, A. Kicheva, C. B¨okel,M. Gonz´alez-Gait´an, and F. J¨ulicher, Development ,1137 (2008).[10] J. O. Dubuis, G. Tkaik, E. F. Wieschaus, T. Gregor,and W. Bialek, Proceedings of the National Academy ofSciences , 16301 (2013).[11] J. Buceta, Journal of The Royal Society Interface (2017), 10.1098/rsif.2017.0316.[12] H. Berg and E. Purcell, Biophysical Journal , 193(1977).[13] N. van Kampen, “Stochastic processes in physics andchemistry,” (Elsevier, Amsterdam ; London, 2007) Chap.XIV, 3rd ed.[14] N. van Kampen, in AIP Conference Proceedings , Vol. 27(AIP, 1976) pp. 153–186.[15] C. Gardiner, “Handbook of stochastic methods forphysics, chemistry and the natural sciences,” (Springer-Verlag, Berlin Heidelberg, 2004) Chap. 8, 3rd ed.[16] Chapter 2 in J. Garcia-Ojalvo and J. Sancho,
Noise inSpatially Extended Systems (Springer-Verlag, New York,1999).[17] Chapter 8 in P. Kotelenez,
Stochastic Ordinary andStochastic Partial Differential Equations: Transitionfrom Microscopic to Macroscopic Equations (Springer-Verlag, New York, 2008).[18] Chapter 1 in H. Holden, B. Øksendal, J. Ube, and T. Zhang,
Stochastic Partial Differential Equations: AModeling, White Noise Functional Approach (Springer-Verlag, New York, 2010).[19] Sec. 1.2.5 in S. V. Lototsky and B. L. Rozovsky,
Stochas-tic Partial Differential Equations (Springer InternationalPublishing, Gewerbestrasse 11, 6330 Cham, Switzerland,2017).[20] Sec. I.1 in L. Landau and E. Lifshitz,
Fluid Mechanics ,Course of Theoretical Physics, Vol. 6 (Elsevier Science,Cambridge, 2014).[21] Sec. 1.2, p.6 in M. Kardar,
Statistical Physics of Fields (Cambridge University Press, Cambridge, 2007).[22] B. Mortazavi, O. Benzerara, H. Meyer, J. Bardon, andS. Ahzi, Carbon , 356 (2013).[23] M. Kojic, M. Milosevic, N. Kojic, K. Kim, M. Ferrari,and A. Ziemys, Computer methods in applied mechanicsand engineering , 123 (2014).[24] D. J. Evans and G. P. Morriss, Statistical Mechanics ofNonequilibrium Liquids (ANUE Press, Canberra, 2007).[25] O. Grimm, M. Coppey, and E. Wieschaus, Development , 2253 (2010).[26] G. Stefanou, Computer Methods in Applied Mechanicsand Engineering , 1031 (2009).[27] K. Sepahvand, S. Marburg, and H. Hardtke, Interna-tional Journal of Applied Mechanics , 305 (2010).[28] Chapter 2 in R. Ghanem and P. D. Spanos, Stochastic Fi-nite Elements: a spectral approach (Dover Publications,US, 31 East 2nd Street, Mineola, NY, 11505, 2003).[29] Sec. I.2 in L. Landau and E. Lifshitz,
Statistical physicsPart 1 , 3rd ed., edited by E. Lifshitz and L. Pitaevskii,Course of Theoretical Physics, Vol. 5 (Pergamon Press,Oxford, 1989).[30] A. Jacobo and A. Hudspeth, Proceedings of the NationalAcademy of Sciences , 15444 (2014).[31] Y. Dublanche, K. Michalodimitrakis, N. K¨ummerer,M. Foglierini, and L. Serrano, Molecular Systems Bi-ology , 41 (2006).[32] C. H. Edwards and D. E. Penney, Elementary Differen-tial Equations , 6th ed. (Pearson Education, New Jersey,2008).[33] Chapter I in S. Chandrasekhar, Rev. Mod. Phys. , 1(1943); J. Piasecki, Acta Physica Polonica Series B ,1623 (2007). [34] M. Kardar, Statistical Physics of Particles (CambridgeUniversity Press, Cambridge, 2007).[35] G. B. Arfken, H. J. Weber, and F. E. Harris,
Mathe-matical Methods for Physicists , 7th ed. (Academic Press,Oxford, 2013).[36] Secs. 2.4 and 4.2 in D. G. Duffy,
Green’s functions withapplications , 1st ed., Studies in advanced mathematics(Chapman & Hall/CRC, 2000 N.W. Corporate Blvd.,Boca Raton, Florida 33431., 2001); Secs. 3.4 and 5.2in
Green’s Functions with Applications , 2nd ed., AppliedMathematics (Chapman and Hall/CRC, 2000 N.W. Cor-porate Blvd., Boca Raton, Florida 33431., 2015).[37] D. Chandler,
Introduction to modern statistical mechan-ics (Oxford University Press, Oxford, 1987).[38] R. Belousov and E. G. D. Cohen, Phys. Rev. E , 062124(2016). [39] L. N. Trefethen, Spectral methods in MATLAB (SIAM,Philadelphia, 2000).[40] F. van de Vosse and P. Minev,
Spectral element methods :theory and applications , Tech. Rep. (Eindhoven Univer-sity of Technology, Eindhoven, 1996).[41] J. P. Boyd,