Diffusive wave dynamics beyond the continuum limit
DDiffusive wave dynamics beyond the continuum limit
Paul B. Dieterle and Ariel Amir ∗ Department of Physics, Harvard University, Cambridge, MA 02138, USA John A. Paulson School of Engineering and Applied Sciences,Harvard University, Cambridge, MA 02138, USA
Scientists have observed and studied diffusive waves in contexts as disparate as population geneticsand cell signaling. Often, these waves are propagated by discrete entities or agents, such as individualcells in the case of cell signaling. For a broad class of diffusive waves, we characterize the transitionbetween the collective propagation of diffusive waves – in which the wave speed is well-describedby continuum theory – and the propagation of diffusive waves by individual agents. We show thatthis transition depends heavily on the dimensionality of the system in which the wave propagatesand that disordered systems yield dynamics largely consistent with lattice systems. In some systemdimensionalities, the intuition that closely packed sources more accurately mimic a continuum canbe grossly violated.
Partial differential equations (PDEs) pervade quanti-tative studies of biological phenomena. One famous suc-cess story of PDEs in biology is the study of diffusivewaves, first begun by Luther in 1906 [1] and furthered byFisher and Kolmogorov et al. in 1937 [2, 3]. In the typi-cal setting, a homogeneous environment can propagate adiffusive wave by picking up some diffusing element – anorganism or a signaling molecule – and producing more ofthat same element – through reproduction or secretion ofthe same entity. The process of diffusion and reproduc-tion can be represented by a partial differential equation,the dynamics solved by way of a traveling wave ansatz,and the wave speed deduced analytically (in certain sce-narios). This type of machinery has been adopted tounderstand traveling waves governing apoptosis [4], eco-logical range expansions [2, 3, 5–7], cell signaling [8–10],mitosis [11], and microbial consortia [12, 13].However, in many of these applications, one is led tobelieve that the underlying PDE models hold in the caseof discrete sources. For instance, cell signaling waves arepropagated by discrete cells [8, 9], mitotic and calciumwaves in Xenopus are paced by discrete sources (nucleiand receptor channels, respectively) [10, 14–18], and re-production in mobile populations happens in centralizedlocations like cities. In such cases, it is not necessarilyvalid to assume a homogeneous environment when cal-culating the properties of the diffusive wave; one mustinstead solve the dynamics for a different set of PDEs –ones with discrete sources.Here, for a broad class of models, we show that thediscreteness of source terms can drastically alter the dif-fusive wave dynamics. We find that in many cases, onecan characterize these corrections analytically for bothordered and disordered systems. The corrections can bedramatic – especially in the case of Fisher wave dynam-ics. In the extreme limit of very disparate sources, wefind the Fisher wave speed is independent of the sourcedensity and the diffusion constant. By examining the ef- ∗ [email protected] fects of discreteness, we see a clear distinction betweendiffusive waves that can be driven by activation at or be-yond the wave front (termed “pulled waves” in the litera-ture) and those that rely on concentration build up fromsources behind the front (termed “pushed waves”). Inthe extreme limit of sparse sources, we show that pulledwaves are driven by self-amplification while pushed wavesrely on excitation from a nearest neighbor. We find thatthe effects of discreteness on asymptotic wave propaga-tion are dimension-dependent. For systems with diffusionin two dimensions, packing sources ever closer does notmake the system behave more like a continuum; for dif-fusion in three dimensions, packing sources ever closermakes the system behave less like a continuum. We also,for the first time, characterize the effects of disorder ondiffusive wave propagation, finding that for many systemdimensionalities a simple lattice model suffices for under-standing the dynamics of diffusive waves in disorderedsystems. MODEL CONSTRUCTION
To begin, we consider a model of diffusive wave propa-gation in which we monitor the concentration, c , of someelements – e.g., signaling molecules in the case of cellsignaling, organisms in the case of Fisher waves – whichdiffuse with diffusion constant D . In addition to basic dif-fusion, the elements can be secreted by the backgroundaccording to some concentration-dependent productionfunction F ( c ). Combining production with diffusion, wearrive at a generic model for the propagation of the con-centration in time: ∂c∂t = D ∇ c + F ( c ) . (1)Concretely, one could imagine a background of cellsof density ρ which secrete diffusible molecules at somerate, a , above a certain threshold concentration, C th ; inthis case, with Θ[ . ] the Heaviside step function, F ( c ) = aρ Θ[ c − C th ] [9, 19, 20]. Inspired by ecology, Fisher a r X i v : . [ n li n . PS ] J a n and Kolmogorov et al. considered production accord-ing to a logistic growth process in which F ( c ) = ac (1 − c )[2, 3]. Importantly, the dynamics driven by the Fisher-Kolmogorov production function depend only on the factthat F ( c ) ∼ c for c (cid:28)
1. We will consider a class of func-tions that interpolate between the Fisher-Kolmogorovand Heaviside limits, in which the production functionis of a Hill function form, F ( c ) = aρ c n c n + C n th . For n = 1,the production function scales linearly in c for c (cid:28) C th : F ( c (cid:28) C th ) ∼ c . In this limit, as we will show shortly,our model reproduces Fisher-Kolmogorov dynamics. For n → ∞ , F ( c ) = aρ Θ[ c − C th ] and we will reproduce Heav-iside dynamics. In total, then, we commit ourselves tostudying equations of the form: ∂c∂t = D ∇ c + aρ c n c n + C n th (2)along with the discrete analogs thereof.Several prior works inspired by calcium waves [15–18]have studied the one dimensional discrete analogs of eq.(1), accounting additionally for simple decay of c by wayof a term − γc . These works cover a variety of produc-tion functions, primarily focusing on the form put forthin eq. (2) as n → ∞ . One such work treats discretenessas a perturbation to the continuum theory dynamics [15].We will instead opt for a solution that yields the wavespeeds – albeit through transcendental equations – evenfar from the continuum limit. (Additional discrete cor-rections have famously been studied by Brunet and Der-rida, among others, in the context of Fisher Waves [21].)Here, we will assume that decay is negligible and willseek to understand the corrections due to discretenessfor all n . Moreover, whereas these works focused almostexclusively on the dynamics of lattice-like systems in onedimension, we will additionally characterize the dynam-ics under disorder and in all dimensions.To take into account the effects of discreteness, we willassume that sources are point-like and thus that the dis-crete analogs of eq. (2) can be written using Dirac deltafunctions δ ( . ). Thusly, and by replacing the density ρ with a sum over the locations r j of all sources, we have ∂c∂t = D ∇ c + a c n c n + C n th (cid:88) j δ ( r − r j ) . (3)The goal of this work will be to study the relationshipbetween the dynamics of eq. (2) and those of eq. (3).To start, we will review the standard, continuum theorywave dynamics of eq. (2) in all dimensions; we will thencalculate the one-dimensional discrete model correctionsto the continuum theory by considering eq. (3) for all n and with different distributions of r j ; finally, specializingto the case of n → ∞ , we will solve for the dynamicsin arbitrary dimensions, showing in the process that thenature of the wave speed corrections heavily depend onthe dimensionality of the diffusive environment and thedistribution of sources within it. CONTINUUM MODEL SOLUTIONS IN ALLDIMENSIONS
Let us begin our study of the continuum model of eq.(2) by considering a purely one-dimensional system, inwhich case we have ∂c∂t = D ∂ c∂x + aρ c n c n + C n th (4)with ρ a density of cells in one dimension. By making atraveling wave ansatz c ( x, t ) = c (˜ x = x − vt ), in whicha concentration wave travels with constant speed v and˜ x is the distance with respect to the wave front, one caneliminate the time derivative of eq. (4) and work insteadwith non-linear ordinary differential equations (ODEs) ofthe form 0 = v ∂c∂ ˜ x + D ∂ c∂ ˜ x + aρ c n c n + C n th . (5)Our goal is to find the wave speed v as a function of thephysiological parameters n , D , C th , a , and ρ .As n → ∞ , this task has a well-known solution.For the latter, one may substitute c n / ( c n + C n th ) → Θ[ c − C th ]; recognize that ˜ x = 0 is the point at which c (˜ x ) = C th ; and stitch together solutions for ˜ x > x <
0, inwhich Θ[ c − C th ] = 1 and 0, respectively [9, 16, 19, 20].Throughout this paper, we will refer to the wave speedwith sources distributed in N dimensions, with diffusionin M dimensions, and with Hill coefficient n as v N,M,n ;the calculation described above reveals that v , ,n →∞ = (cid:112) aρD/C th . (6)The solution for n = 1 is more involved and appeals tomathematics first worked out by Fisher and Kolmogorov et al. [2, 3]; here, one examines the region well be-yond the wave front (˜ x (cid:29) D/v ), in which c (cid:28) C th and c/ ( c + C th ) ≈ c/C th . In these limits, we may consider amodified form of eq. (5):0 ≈ v ∂c∂ ˜ x + D ∂ c∂ ˜ x + aρc/C th . (7)One then plugs an ansatz of the form c = c e − γ ˜ x into eq.(7), which results in a relationship between v and γ : v = Dγ + aρ/C th γ. (8)Given the restriction that v, γ >
0, eq. (8) appears to al-low a range of solutions v ≥ (cid:112) aρD/C th ; famously [2, 3],the minimum speed allowed by eq. (8) is the one that thewave attains, a fact that has been proven to hold with pe-riodic and homogeneous production functions alike [22].Thus, v , , = 2 (cid:112) aρD/C th . (9)It is not a coincidence that v , , and v , ,n →∞ bothscale as (cid:112) aρD/C th . As has been shown previously[1, 19], this scaling can be deduced by dimensional anal-ysis considerations. In brief, the source terms on theR.H.S. of eq. (4) are proportional to aρ , meaning allconcentration scales are also proportional to this factor.Thusly, the four parameters of a , ρ , C th , and D are boileddown to two parameters, D and C th /aρ with respectivedimensions m /s and s. As we wish to calculate an emer-gent wave speed v with units m/s, the only viable formulafor v is a scaling law of the form v ∼ (cid:112) aρD/C th .While these results plainly hold in systems of dimen-sionality ( N, M ) = (1 , N = M , so long as thedynamics are observed at radii r (cid:29) N D/v [6, 19].In contrast to systems with N = M , one could con-sider a system in which some N -dimensional distribu-tion of sources is embedded within a higher-dimensional, M > N diffusive environment. One example is an assem-blage of cells sitting upon a two-dimensional plane with asemi-infinite extracellular medium above [19]. In such ascenario, assuming a semi-infinite diffusive environmentwith N = 1 and M = 2, eq. (2) becomes ∂c∂t = D (cid:18) ∂ c∂x + ∂ c∂z (cid:19) + 2 aρδ ( z ) c n c n + C n th (10)with x the dimension in which the sources live, z theperpendicular dimension, and ρ a one-dimensional den-sity of sources. Within the aforementioned traveling waveansatz, one can reduce this to a spatial PDE of the form0 = v ∂c∂ ˜ x + D (cid:18) ∂ c∂ ˜ x + ∂ c∂z (cid:19) + 2 aρδ ( z ) c n c n + C n th . (11)Before noting the exact solutions of the wave speed, weagain appeal to dimensional analysis. As in the case of N = M discussed above, the source terms on the R.H.S.of eq. (4) are proportional to aρ and thus we are again leftwith only two effective model parameters, D and C th /aρ ,with which to construct a wave speed. However, unlikethe case of N = M , in which C th /aρ ∼ s, we now have C th /aρ ∼ s/m. Thus, the only formula involving D and C th /aρ which forms a valid wave speed, v , of units m/sis the D -independent scaling law of v ∼ aρ/C th .For n → ∞ , one can solve eq. (11) using a partialFourier transform over z and the methodology employedon eq. (5) to yield v , ,n →∞ = 2 aρ/πC th . (12) N = M N = M − n = 1 v N,N, = 2 (cid:112) aρD/C th v N,N +1 , = 2 aρ/C th n → ∞ v N,N,n →∞ = (cid:112) aρD/C th v N,N +1 ,n →∞ = 2 aρ/πC th TABLE I.
Summary of continuum theory wave speeds:
Here, we summarize the diffusive wave speed for systems withsources in N dimensions and diffusion in M dimensions. Weconsider system dimensionalities of N = M and N = M − n = 1 and n → ∞ . As shown in the main text,dimensional analysis demands that systems with N = M − D ) independent wave speed andscale differently in the source density ρ , emission rate a , andthreshold concentration C th . A similar approach with n = 1 requires an unwieldy par-tial Fourier transform over z of0 = v ∂c∂ ˜ x + D (cid:18) ∂ c∂ ˜ x + ∂ c∂z (cid:19) + 2 aρδ ( z ) c (˜ x, z ) /C th . (13)In Appendix D we show it is possible to solve for thedynamics of eq. (13) with Green’s function methods.Doing so, we show that v , , = 2 aρ/C th . (14)In closing this section we note that, as in the case of N = M , eq. (12) and eq. (14) hold for all systemswith N = M − r (cid:29) N D/v ). We also note that systems forwhich N = M − j and j ≥ v increases ( j >
2) orremains constant ( j = 2) with increasing C th . The is-sue is that embedding a continuum of diffusive sourcesin a higher-dimensional space creates a divergent con-centration at the sources when the discrepancy betweenthe source dimension and the diffusive dimension is morethan one.The results for continuum systems with N = M or N = M − DISCRETE MODEL SOLUTIONS IN ONEDIMENSION
We now seek to depart from the continuum theory ofdiffusive waves by focusing on waves propagated by dis-crete sources. To do so, we examine the one-dimensionalversion of the discrete model put forth in eq. (3): ∂c∂t = D ∂ c∂x + a c n c n + C n th (cid:88) j δ ( x − x j ) . (15)Here, we have discrete sources at positions x j , each ofwhich responds to the local concentration c by emittingdiffusible molecules at some rate ac n / ( c n + C n th ). In thissetting, we would like to find something akin to the wavespeed calculated in the continuum limit.To begin, let us consider the dynamics of sources on alattice with separation d (i.e., density ρ = 1 /d ) as in Fig.1A/B. We anticipate that the wave speed will be closeto that of the continuum model for sufficiently small d –but how small?Examining the n = 1 solution above, we see that ouransatz has a beyond-the-front solution of e − vx/ D . Thus, D/v is roughly the length scale of the diffusive wave,and we expect that when many sources fit inside thislength scale – i.e., when dv/D (cid:28) v itself dependsupon d . Within the continuum theory ρ = 1 /d = ⇒ v ∼ (cid:112) aD/dC th , which implies that d (cid:28) DC th /a is the limitin which the two models should agree and that DC th /a is the length scale against which we should compare d .We can come about this length scale through simplermeans of dimensional analysis. Omitting d , the onlyother parameters in our model are C th , a , and D . Asin the continuum limit, we may eliminate a by notingthat the source terms in eq. (15) are all proportional to a . Thus, we are left with D (units of m /s) and C th /a (units of s/m ) as the only other model parameters. Usingthese parameters to construct a length scale, we see thereis a unique length scale DC th /a against which we maycompare d . In essence, when C th /a or D is sufficientlysmall, only the nearest neighbors contribute to the lo-cal concentration and we can expect deviations from thecontinuum theory. Keener colorfully refers to this limitas the “saltatory” regime, in which the wave front hopsform one source to the next [15]. With this length scalein mind, let us conduct a rigorous analysis for n → ∞ and n = 1. Lattice solution for n → ∞ In the limit n → ∞ , it is possible to construct a rela-tionship which gives a unique wave speed v as a functionof d , D and C th /a . To do so, we consider a lattice ofpoint sources, as in Fig. 1B, through which a diffusivewave has been traveling at speed v . We set the time andspace origins such that the concentration at t, ˜ x = 0 is C th ; doing so, we stipulate that the source a distance jd away has been continuously emitting at a rate a sincetime t = − jd/v . Using the diffusion equation Green’sfunction for a point-like source, we can therefore see thatthe concentration, c j , created by the source a distance jd away from the origin is given by c j = a (cid:90) − jd/v d ˜ t e j d / D ˜ t / (cid:112) − πD ˜ t = jad D (cid:32) e − jdv/ D (cid:115) Dπjdv + erf (cid:114) jdv D − (cid:33) (16)while the threshold concentration, C th , obeys a self-consistency relationship given by summing over all thesources: C th = ∞ (cid:88) j =1 c j . (17)In constructing eq. (17), we have successfully derived arelationship between v and the system parameters C th /a , d , and D . We will now explore this relationship to testthe intuition we developed through dimensional analysis.By examination of eq. (16), we plainly see that thelength scale with which to compare d is D/v . The latteris the length scale of the diffusive wave, as mentionedabove. Let us assume that d (cid:28) D/v – i.e., that manysources fit within the length scale of the diffusive wave –in which case one may turn eq. (17)’s sum over j into anintegral and deduce that d (cid:28) D/v : C th ≈ (cid:90) ∞ dj c j = aD/v d = ⇒ v = (cid:112) aD/dC th , (18)which is exactly the continuum theory result of eq. (6)with ρ = 1 /d . Using our expression for v , we may deducethat the assumption of d (cid:28) D/v is equivalent to: d (cid:28) D/v = D (cid:112) dC th /aD = ⇒ d (cid:28) C th D/a, (19)exactly as predicted by our dimensional analysis. In thislimit, the concentration wave can be understood as a col-lective phenomenon in which many sources contribute tothe concentration at the wave front (Fig. 1F) and thewave speed is approximately that of our continuum the-ory (Fig. 1C). Thus, when the dimensionless parameter φ = ad/DC th is small, the collectivity is high and viceversa .In the contrasting limit of d (cid:29) C th D/a , the concen-tration at the wave front is primarily generated by thenearest source (Fig. 1F) and the wave speed is consid-erably less than that predicted by the continuum theory.We will return to this limit shortly.
Lattice solution for n = 1 and numerics for all n Next, we consider the n = 1 limit of eq. (15). Toaddress this limit, we again consider a lattice of point a A B n =1 n !1 c ( x =0, t =0) = C th j =123...emitting since t = - kd / vk vdd N N +1 N -1 ... N + jc ( x =( N + j ) d , t =0) » e -( N + j ) °d v ( ° ) v / v , , n v C t h / a ad / DC th = Á C D E c j / c j , m a x n =1-10 100 F -10 100 n !1 j j ad / DC th = Á
10 2 collectivity n =1 Á =1301000 Á =1301000 n =11.5235... FIG. 1.
Discrete lattice wave dynamics in 1D. A:
Schematic demonstrating our solution method for Fisher-like waveswith n = 1. Here, we examine the region well beyond the wave front and make an ansatz that the concentration at source N + j is given by c ( x = ( N + j ) d, t = 0) ∼ e − ( N + j ) γd . This ansatz allows a range of solutions v ( γ ), the minimum of whichis the realized wave speed. B: Schematic demonstrating the solution method for Heaviside waves with n → ∞ . Here, theconcentration at the wave front is given by C th while the source a distance kd from the wave front has been emitting since time t = − kd/v . C: Discrete wave speeds divided by continuum wave speed as a function of ad/DC th = φ for the n = 1 theory(blue line), n → ∞ theory (black line), and finite n numerics (colored circles). For all n >
1, large separation between sourcesleads to a sub-continuum wave speed while for n = 1 the wave speed increases compared with its continuum value. At smallvalues of φ , the wave can be understood as a collective effect that is well-described by a continuum theory. D: Theoretical wavespeed for n = 1 as stipulated by eq. (23). At large values of φ , the wave speed attains a d - and D -independent value of a/ C th . E: Concentration at source N generated by source N + j for various values of φ and for n = 1. At small φ , many sourcescontribute to the concentration at source N while as φ increases, the wave propagation becomes a few-body phenomenon. For φ = 1000, we see the extreme limit of a pulled wave, in which the concentration at the wave front is almost entirely due toself-amplification. F: Concentration at the wave front generated by source j as n → ∞ . As with the case of n = 1, many sources contribute tothe concentration at the wave front when φ is small. As φ increases, the wave propagation becomes a few-body problem. Incontrast to the case of n = 1, for φ = 1000 the wave front is pushed by the nearest neighbor behind the front. sources, as in Fig. 1A, in which a diffusive wave withspeed v has been traveling. Examining the region ofspace well beyond the wave front in which c (cid:28) C th , wefollow the example of Fisher [2] and propose an ansatzthat the concentration some distance ˜ x = N d beyond thewave front is given by: c (˜ x = N d, t ) = c e − γNd + γvt . (20)In the limit of c (cid:28) C th , the production function may beapproximated as c/ ( c + C th ) ≈ c/C th . Thus, the con-centration, c j , generated by the source ˜ x = ( N + j ) d asmeasured at ˜ x = N d and t = 0 can again be calculatedby integrating the production function of source N + j against the Green’s function for the diffusion equation: c j = ac e − γNd C th (cid:90) −∞ dt e j d / Dt e − γjd + γvt √− πDt = ac e − γNd C th √ γDv e − d (cid:16) jγ + | j | √ γv/D (cid:17) . (21)One may construct a self-consistency relationship alongthe lines of eq. (17) by recognizing that the concentration c (˜ x = N d, t = 0) is merely the sum of all the c j , so that c (˜ x = N d, t = 0) = ∞ (cid:88) j = −∞ c j , (22)which can be written in full as2 C th √ Dγva = ∞ (cid:88) j = −∞ exp (cid:104) − jγd − | j | d (cid:112) γv/D (cid:105) =1 + 1 e dγ + d √ γv/D − e − dγ + d √ γv/D − . (23)Through eq. (23), we are endowed with a relationship,like eq. (17), that fully specifies v given the physiologicalparameters C th /a , d , and D – but with one exception.As in the famous case of Fisher’s equation, we must alsospecify γ ; and, as in the case of Fisher’s equation, thevalue of γ that is selected is the one that minimizes v .Given eq. (23), one can find the wave speed v by min-imizing v over all γ subject to the transcendental rela-tionship of eq. (23). That the wave speed achieved isthe minimum has been proven for all periodic systems inbeautiful recent work [22].As in the case of n → ∞ , our goal here is to see howthe wave speed stipulated through the above calculationsrelates to the continuum theory of eq. (9). Appealing toour previous intuition, we can guess that d (cid:28) C th D/a will yield a wave speed consistent with the continuumtheory. Indeed, we may expand eq. (23) to lowest orderin d for d (cid:28) γd, d (cid:112) γv/D ; doing so leaves one withexactly the relationship of eq. (8) with ρ = 1 /d . Thevalue of γ which minimizes v is γ = (cid:112) dC th D/a , meaningthat our imposition of d (cid:28) γd, d (cid:112) γv/D is equivalent to d (cid:28) C th D/a . In this limit, the wave speed v agrees wellwith the continuum theory (Fig. 1C) and the effect isa collective one in which many sources contribute to theconcentration generated at every other source (Fig. 1E).In the opposite limit of d (cid:29) C th D/a , the primarydriver of the concentration generated at each source is theconcentration generated by the source itself (Fig. 1E).The source is seeded with some small concentration fromits nearest neighbor, then self-amplifies; thus, unlike inthe case of n → ∞ in which only the nearest neighbormatters, we now must account for both the nearest neigh-bor and the self-interaction. Surprisingly, the wave speedin this limit is independent of d and D , as we will nowshow and as is pictured in Fig. 1D. In this limit, one canapproximate the sum in eq. (23) by considering only thedominant j = 0 , − d (cid:29) C th D/a : 2 C th √ Dγva = 1 + e γd − d √ γv/D . (24)Treating v as a function v ( γ ), one can take a derivativeof both sides of the above equation with respect to γ .Rearranging terms, we can see that ∂v/∂γ = 0 in the d → ∞ limit requires γ − (cid:112) γv/D ∼ log d/d and thusthat d (cid:29) C th D/a : γ = v/D = ⇒ v = a/ C th ,γ = a/ DC th (25) where the relationship after the arrow is had by pluggingthe relationship before the arrow into eq. (24).It is strange that, in the d (cid:29) C th D/a limit, a wavepropagated between discrete sources and driven by dif-fusion has a speed that depends neither on the distancebetween the sources nor on the rate of diffusion. To un-derstand how this can be, we return to the ansatz of eq.(20) and plug in the values of eq. (25), which yields d (cid:29) C th D/a : c (˜ x = N d, t ) = c e − Nda DC th + a t DC . (26)An examination of eq. (26) shows that the wave speedis indeed d - and D -independent in the large- d limit. Wesee that c (˜ x = N d, t ) = c [˜ x = ( N + 1) d, t + 2 C th d/a ],meaning that the concentration wave travels a distance d in a time τ = 2 C th d/a . As a result, the wave speed v = d/τ is independent of d . In essence, increasing d resultsin an exponentially smaller initial concentration at anypoint source but gives that point source more time to self-amplify – effects that exactly cancel when calculating thewave speed.Similarly, increasing the diffusion constant, D , resultsin an exponential increase in the initial concentration ata point source (see the factor of e − Nda/ DC th in eq. (26))due to the fact that molecules can more quickly diffusefrom source to source; however, this effect is exactly can-celled by the decrease in the concentration’s growth rate, a / DC , and results in the D -independent hoppingtime between sources, τ = 2 C th d/a , calculated above.To understand the origin of this d - and D -independence, we must therefore understand why theconcentration takes the precise form of eq. (26). Todo so, let us consider a single source seeded with somesmall concentration c (cid:28) C th at t = 0 and secretingwith a concentration-dependent rate ac/C th . In this sce-nario, the concentration, c , at the source at time t is givenself-consistently given by integrating ac/C th against thediffusion equation Green’s function, thus giving us c ( t ) = aC th √ πD (cid:90) t d ˜ t c (˜ t ) / (cid:112) t − ˜ t, (27)which can be solved with an ansatz c e Γ t for t (cid:29) / Γ bysending t → ∞ in the integral bound, telling us that theconcentration at the source indeed grows exponentiallywith rate Γ = a / DC . (28)Comparing this exponential growth rate of the concen-tration at a single point source with the traveling wavedynamics stipulated by eq. (26), we can plainly see thatthe growth rate of the concentration at sources well be-yond the wave front is exactly the exponential growthrate of the concentration at a single, self-amplifying pointsource.Moreover, by integrating the Green’s function againstthe concentration production rate of ac/C th for c (cid:28) C th ,we can deduce that a self-amplifying point source thathas been growing for a time t (cid:29) Γ , r /D creates a con-centration profile given by c ( x, t ) = aC th √ πD (cid:90) t d ˜ t e − x / D ( t − ˜ t ) e Γ˜ t (cid:112) t − ˜ t ∼ (cid:90) ∞ d ˜ t e − x / D ( t − ˜ t ) e Γ˜ t (cid:112) t − ˜ t ∼ e Γ t − x √ Γ /D . (29)This is precisely the profile that the diffusive wave gener-ates in the limit d (cid:29) C th D/a , see eq. (26). Thus, we canunderstand the temporal dependence of eq. (26) as theself-amplification rate of a single source and the spatialdependence of eq. (26) as the the concentration profilethat results from this self-amplification.
Transition to pushed waves for n > While in the case of n = 1 we may solve for the dy-namics through a Green’s function analysis of eq. (15) inthe c (cid:28) C th limit, the same analysis does not permit so-lutions for n > n = 1 . n > n = 1, as is well-known fromthe studies of continuum models like eq. (4), in which n = 1 waves are termed “pulled” by self-amplification inthe c (cid:28) C th limit while n > n > c (cid:28) C th . In this limit, c n / ( c n + C n th ) ≈ c n /C n th and we may construct a self-consistency relationshipof the form eq. (27) by integrating the emission rate ac n /C n th against the Green’s function for the diffusionequation and setting that equal to the local concentra-tion at some later time t : c ( t ) = aC n th √ πD (cid:90) t −∞ d ˜ t c n (˜ t ) / (cid:112) t − ˜ t, (30)Fixing the time origin such that c ( t = 0) = c , theconcentration grows as (with α ∼ − a / C n th D , see Ap-pendix A for full expression) c ( t ) = (cid:104) αt + c − n − (cid:105) − n − , (31)meaning that, for c (cid:28) C th , the time for c ( t ) to doublescales as t ∼ ( C th /c ) n − . (32) In contrast, for any n , a diffusive wave (whether propa-gated by discrete or continuous sources, see Appendix B)will create a concentration profile in the limit c (cid:28) C th that decays exponentially in distance beyond the wavefront and thus grows exponentially in time. Therefore,the time scale needed for the wave to pass through a dis-tant point source with initial concentration c is roughlygiven by: t ∼ log [ C th /c ] . (33)We can therefore see that for n >
1, individual point-sources cannot propagate an initial concentration quicklyenough to exceed C th by the time the wave passesthrough; thus, the concentration at sources well beyondthe wave front (where c (cid:28) C th ) is primarily generatednot through the self-interaction of sources in that re-gion but rather through the build up of concentrationgenerated by sources behind or around the wave front.This is precisely the phenomenon of a “pushed” wave.This distinction is shown clearly in Figs. 1E/F. Forlarge ad/DC th , we can see that the concentration at thewave front for n = 1 is generated almost entirely by self-amplification. In contrast, the concentration at the wavefront for n → ∞ is generated almost solely by the sourceimmediately behind the wave front.Because the wave front for systems with n > n → ∞ in that the wave front will hop to the sourcebeyond the wave front primarily as a result of concentra-tion build up generated by the source immediately insidethe wave front. This is a necessity because sources with n > τ ∼ d to bring a neighboring source over the threshold concen-tration. The apparent wave speed of this process is then v ∼ d/τ ∼ /d . In comparison, the continuum theorywave speed scales as v , ,n ∼ / √ d . Thus, for n > ad/DC th (cid:29)
1, the discrete source wave speed will fall be-low the continuum theory wave speed as ad/DC th → ∞ .To affirm this result, we turned to numerical simula-tion. In Fig. 1C, we plot the wave speeds for various val-ues of n and ad/DC th . We obtained these values throughnumerical simulation of eq. (15) and compared them totheir continuum theory analog through numerically solv-ing eq. (4) (see Appendix C for numerical details). Asexpected, systems for all values of n agree well with theircontinuum theory analogs for small ad/DC th . For large ad/DC th , systems with n = 1 display a super-continuumwave speed; in contrast, for sufficiently large ad/DC th ,systems with n > A v = d / ¿j =1 2 3 4 ... N -1 N N +1 ¿ ¿ ¿ ¿ N -1, N ¿ N , N +1 v / v , , n n =1 n !1 -1 B C ad / DC th ad / DC th FIG. 2.
Comparison of lattice and disordered wavespeeds in 1D. A:
Solution methodology for finding thewave speed propagated by a discrete, disordered system in1D. The wave speed is determine by the average distance be-tween sources, d , divided by the mean hopping time, ¯ τ . In thelimit d (cid:28) DC th /a , ¯ τ can be solved for analytically. B: Dis-crete wave speed for n = 1 divided by the continuum theorywave speed as a function of ad/DC th . When v/v , ,n ≈ ad/DC th (cid:29) C: Same as B , but for n → ∞ .Unlike for n = 1, here we see significant deviation of thedisordered numerics (circular markers) and the lattice theory(solid line). In the limit ad/DC th (cid:29)
1, this discrepancy iswell-characterized by our nearest neighbor disordered theory(dashed grey line). as they rely not on self-excitation but on the compara-tively slower process of diffusion between discrete sourcesaround the front.
Effects of disorder for n = 1 , n → ∞ Of course, in biology and elsewhere, the sources whichpropagate diffusive waves are not always found in discretelattices. As such, we now seek to quantify how similarthe dynamics of discrete lattice systems are to the dy-namics of disordered relays with discrete sources. Here,we will specialize to the cases of n = 1 and n → ∞ asthese limits will afford us analytical expressions for thewave speed of disordered systems for ad/DC th (cid:29)
1. Asin the case of lattice systems, the relevant dimensionless parameter for disordered systems is ad/DC th ; this factcan be gleaned through dimensional analysis or by con-sidering the concentration variance generated by a fixedspeed diffusive wave propagating through a disorderedmedium (for details of the latter, see Appendix E).In a disordered system, a chain of sources in 1D willhave some inhomogeneous density of sources (Fig. 2A).The result of this straightforward observation is that thewave front does not propagate smoothly through the sys-tem of discrete sources, but rather hops incongruouslyfrom source j to source j + 1 in a time τ j,j +1 (Fig. 2A).One can therefore define the wave speed, v , in a disor-dered 1D system only through an ensemble average inwhich the average number of sources with the local con-centration exceeding threshold at time t , (cid:104) n ( t ) (cid:105) , is relatedto the wave speed v via the average source separation, d : (cid:104) n ( t ) (cid:105) = vt/d (34)where we assume sufficiently large t so that the initialcondition is negligible.In general, it is difficult to fully perform the calcula-tion for (cid:104) n ( t ) (cid:105) ; to see why, consider the case of n → ∞ and consult Fig. 2A. In order to find the hopping timebetween source j and source j + 1, one must know theconcentration at source j + 1 at all times so that one canfind the moment at which it exceeds C th . But this re-quires knowledge of when sources j , j −
1, etc. turnedon (i.e., first exceeded C th ), so that calculating τ j,j +1 requires knowing every previous hopping time.For all n , one can break this endless cycle in the limit ad/DC th (cid:29) τ j,j +1 is a function (to be calculated shortly) only of n , D , C th /a , and the distance, d j,j +1 , between source j andsource j + 1. Thus, the hopping time between sourcesseparated by a distance x , τ ( x ) can be written in termsof the lattice wave speed, v lattice1 , ,n : τ ( x ) = x/v lattice1 , ,n . (35)Then, given some distribution of separations, p ( x ), withmean d , we can calculate (cid:104) n ( t ) (cid:105) as t divided by the meanhopping time, ¯ τ : d (cid:29) DC th /a : (cid:104) n ( t ) (cid:105) = t/ ¯ τ = t (cid:18)(cid:90) ∞ dx τ ( x ) p ( x ) (cid:19) − . (36)It follows from eq. (34) that d (cid:29) DC th /a : v = d/ ¯ τ = d (cid:18)(cid:90) ∞ dx τ ( x ) p ( x ) (cid:19) − . (37)For Poisson-distributed sources with mean separation d , p ( x ) = e − x/d /d and it only remains to calculate τ ( x ), atask we shall now undertake for n = 1 and n → ∞ .Our task is straightforward for n = 1 since v = a/ C th when d (cid:29) DC th /a . Thus, τ ( x ) = 2 C th x/a and n = 1 , d (cid:29) DC th /a : v = d/ ¯ τ = a/ C th (38)by eq. (37). Note that this result holds for all p ( x ) as d = (cid:104) x (cid:105) and ¯ τ = C th a (cid:82) dx xp ( x ) = C th (cid:104) x (cid:105) a . This dis-ordered wave speed is plotted (dashed line) in Fig. 2Balong with the lattice theory (solid line) and results fromnumerical simulations of Poisson-disordered systems (cir-cles). As seen in Fig. 2B, the lattice theory closely cor-responds with numerical simulations over several ordersof magnitude in ad/DC th .For n → ∞ , C th can be approximated by truncatingeq. (17) at j = 1. We then refer to eq. (16) and notethat on a lattice v = d/τ , meaning that n → ∞ , d (cid:29) DC th /a : C th = ad D (cid:32) e − d / Dτ (cid:114) Dτπd + erf (cid:114) d Dτ − (cid:33) , (39)a relationship that can be numerically inverted for all d ,then numerically integrated over to give v .The resulting nearest neighbor theory for v (dashedline) is plotted in Fig. 2C along with the lattice the-ory (solid line) and results from numerical simulations ofPoisson-disordered systems (circles). Here, we see sub-stantial deviation of the lattice theory and the numericalsimulations of disordered systems for intermediate andlarge values of ad/DC th . However, the disordered sys-tem numerics agree very well with our nearest neighbortheory in the limit ad/DC th (cid:29) n = 1 and all values of ad/DC th .In contrast, with sufficiently large ad/DC th and n → ∞ ,Poisson-disordered and lattice systems propagate diffu-sive waves of substantially differing wave speed. DISCRETE MODEL SOLUTIONS IN HIGHERDIMENSIONS
We have fully characterized the diffusive wave dynam-ics for sources and diffusion in one dimension ( N = M =1) and now turn our attention to systems of other dimen-sionalities. Interestingly, we will see that whereas, e.g.,the N = M continuum theories all give the same asymp-totic wave dynamics, the discrete N = M = 1 dynamicsare distinct from those of N = M = 2 ,
3; similarly, thedynamics of N = M − N = M − M ≥
2. This is due to the fact thatthe local concentration at a point source diverges intwo-or-more dimensions. This results in an infinite self-interaction. Mathematically, one can see this by inte-grating the Green’s function of the diffusion equation fora continuously emitting point source in M dimensions: M ≥ c ( x = 0 , t ) = a (cid:90) t d ˜ t (cid:2) πD ( t − ˜ t ) (cid:3) − M/ → ∞ . (40)Therefore, the self-interaction will diverge for point-likesources with M ≥ n → ∞ . In this limit, the dynamics are well-posedbecause sources do not interact with their own emissionsuntil they are above the c = C th threshold, at whichpoint self-interaction is irrelevant.We will momentarily perform the task of solving forthe wave speeds v within lattice and Poisson-disorderedarrangements, but let us first develop intuition for thesecases through dimensional analysis. We have committedourselves to studying models of the form ∂c∂t = D ∇ c + a Θ[ c − C th ] (cid:88) j δ ( r − r j ) , (41)in which we recognize (as we have previously) two inde-pendent parameters, D and C th /a , along with the meanseparation between sources d .As above, our task here is to understand how varying d affects the wave speed v for given values of D and C th /a .For M = 1, we have seen that we must compare d to thediffusive wave length scale D/v and that doing so revealsthe continuum theory to be valid when d (cid:28) DC th /a .The same does not hold for M = 2. Here, D/v isstill the length scale of the diffusive wave, but now, with N = 2, v ∼ (cid:112) aD/d C th within the continuum theory.Constructing a dimensionless parameter dv/D (cid:28) d and merely requires that a/DC th (cid:28) d has nothing to do with the agree-ment between the discrete and continuum theories. In-stead, the agreement is entirely determined by the valueof DC th /a . A similar argument holds for ( N, M ) = (1 , d . The parameters we can use – bythe same arguments as in the M = 1 case – are D (units0m /s) and C th /a (units m /s, as distinguished from theunits of m/s for M = 1). Thus, there is no combinationof D and C th /a which can form a length scale. For fixedvalues of D and C th /a , we cannot pack sources closertogether and get better agreement with the continuumtheory. Nonetheless, our previous intuition that D and C th /a ought to be small in order to violate the continuumtheory is still valid. Here, the dimensionless parameterwhich describes the agreement between the discrete andcontinuum theories is a/DC th .The situation is even more extreme for M = 3. Yetagain, we seek to compare d with D/v . From our con-tinuum theory results, we know that when N = 3, v ∼ (cid:112) aD/d C th . Requiring dv/D (cid:28) dDC th /a ) − (cid:28) worse when one decreases d at fixed D and C th /a . Thus, bycramming sources ever closer together in a 3D diffusiveenvironment, one would weaken the agreement betweenthe discrete system dynamics and those of the continuum.A similar argument holds for N = 2. We will shortlydemonstrate this effect through exact calculations andnumerical simulation. Lattice solutions
Our goal in this section is to systematically find therelationships between v , D , d , and C th /a for a set ofpoint sources sitting on a lattice in N dimensions withdiffusion in M dimensions. We have already solved thecase of N = M = 1 above (see eqs. (16) and (17)).Our methodology for the solutions with other N, M islargely the same as that already employed, but requiressome further elaboration.To calculate self-consistency relationships analogous toeq. (17) for N = 2 , N -dimensions.With N = 2, we obtain a picture similar to Fig. 1B,but with sources at y = 0 , ± d, ± d, . . . . At t = 0,the sources at x = − jd and y = 0 , ± d, ± d, . . . havebeen emitting since t = − jd/v . If we assume that c ( x = y = t = 0) = C th , then we can construct aself-consistency relationship by adding up the contribu-tions from all the point sources left of the origin. Forthe sake of clarity, we have performed this calculationfor ( N, M ) = (2 ,
2) below; the self-consistency relation-ship for other system dimensionalities can be found inAppendix F. (N,M)=(2,2):
To find a self-consistency relationshipfor this dimensionality, we consider the concentration, c ( M =2) j,k , generated at x = y = t = 0 by a diffusive pointsource in two dimensions at ( x, y ) = ( − jd, kd ) that hasbeen emitting since t = − jd/v . This can be calculatedby integrating the diffusion equation Green’s function fora 2D environment, a process that yields, with Γ [ . ] theincomplete Gamma function, c ( M =2) j,k = − a πD (cid:90) − jd/v d ˜ t e j d k d D ˜ t / ˜ t =Γ (cid:20) vd Dj ( j + k ) (cid:21) / πD. (42)We may therefore calculate the self-consistency relationfor ( N, M ) = (2 ,
2) by adding up all of the sources:(
N, M ) = (2 ,
2) : C th = ∞ (cid:88) j =1 ∞ (cid:88) k = −∞ c ( M =2) j,k . (43)Using eq. (43) and its analogs [eqs. (69), (71), and(72)], we have plotted the wave speeds, v , for everysystem dimensionality in Fig. 3 (solid lines). As ex-pected, these relationships reveal that for M = 2 thediscrete-continuum correspondence is entirely indepen-dent of d and worsens as a/DC th is increased; for M = 3,the discrete-continuum correspondence worsens as d de-creases. Thus, our previously developed intuition anddimensional analysis provide clear intuition for the dy-namics of diffusive waves propagated by discrete sources. Effects of disorder: theory and numerics
Our goal in this section is to examine the effects ofdisorder on the diffusive wave propagation in higher di-mensional systems. In short, we find that disorder leadsto mild deviations from the lattice wave theory developedin the previous section, with the most severe discrepan-cies for (
N, M ) = (1 ,
2) (Fig. 3).To begin, we note that systems with N ≥ N, M ) = (1 , N, M ) = (1 ,
2) andcan numerically examine the effects of disorder in all sys-tem dimensionalities – tasks we will now undertake.In the non-collective limit of an (
N, M ) = (1 ,
2) sys-tem, we can arrive at an exact expression for the wavespeed by considering the mean hopping time betweennearest neighbors. With v lattice1 , ,n →∞ as the wave speed ofa lattice theory, the hopping time between two sourcesseparated by a distance x is given by: τ ( x ) = x/v lattice1 , ,n →∞ . (44)By our previous dimensional analysis and eq. (12), weknow that v lattice1 , ,n →∞ = f ( a/DC th ) v , ,n →∞ = f ( a/DC th ) 2 aπxC th (45)1 v / v , , n !1 v / v , , n !1 v / v , , n !1 v / v , , n !1 ( N , M )=(1, 2) ( N , M )=(2, 2) ( N , M )=(2, 3) ( N , M )=(3, 3) A C D a / DC th d (a.u.) a / DC th d (a.u.) a / dDC th a / dDC th B FIG. 3.
Comparison of lattice and Poisson-disordered wave speeds in higher dimensions. A:
Wave speed dividedby the continuum theory value with (
N, M ) = (1 ,
2) for disordered systems (circles) and lattice systems (black line); we plotthis quantity for varying a/DC th (top) and varying d (bottom). When v/v N,M,n →∞ ≈
1, a discrete system is in the continuumlimit. As predicted by dimensional analysis, v/v , ,n →∞ does not change at all with d (bottom panel), but decreases as a/DC th increases (top panel). Results in the bottom panel are shown for a/DC th = 10. The disordered and continuum theoriesshow substantial deviation, which can be understood through a disordered nearest neighbor theory (dashed grey line) at large a/DC th . B: Same as A , but for ( N, M ) = (2 , N, M ) = (1 , a/DC th (toppanel) and does not vary with d (bottom panel). Results in the bottom panel are shown for a/DC th = 1000. C: Wave speeddivided by the continuum theory value as a function of a/dDC th with ( N, M ) = (2 , M = 1 , d results in larger deviations from continuum theory. D: Same as C but for ( N, M ) = (3 , a/dDC th . Again, we observe that decreasing d results in larger deviations from the continuum theory, as predicted by dimensional analysis. where f ( a/DC th ) is a function that describes the solidcurve plotted in Fig. 3A. Thus, the average hoppingtime, ¯ τ , is given by¯ τ = (cid:90) ∞ dx p ( x ) τ ( x ) = πC th af ( a/DC th ) (cid:90) ∞ dx x p ( x ) , (46)which for Poisson-distributed sources with mean separa-tion d gives us v = d/ ¯ τ = f ( a/DC th ) aπxC th = v lattice1 , ,n →∞ / . (47)This relationship agrees with numerical simulations ofdisordered systems in the large- a/DC th limit.To understand the effect of disorder in the remainingsystem dimensionalities, we appeal to numerical simula-tion. As shown in Figure 3, the wave speeds observed inPoisson-disordered systems largely agree with the wavespeeds observed in lattice systems of comparable a/DC th ( M = 2) or a/dDC th ( M = 3).As expected from dimensional analysis, the effect ofvarying d in systems with M = 2 is nil (Figs. 3A/B). Thisintuition holds for both Poisson-disordered and latticesystems. Meanwhile, increasing d in systems with M = 3indeed has the counter-intuitive effect of improving theagreement with the continuum theory (Figs. 3C/D). In systems with N >