Dynamics of diffusive cell signaling relays

  1. Paul B Dieterle
  2. Jiseon Min
  3. Daniel Irimia
  4. Ariel Amir  Is a corresponding author
  1. Department of Physics, Harvard University, United States
  2. Department of Molecular and Cellular Biology, Harvard University, United States
  3. BioMEMS Resource Center and Center for Surgery, Innovation and Bioengineering, Department of Surgery, Massachusetts General Hospital, United States
  4. John A. Paulson School of Engineering and Applied Sciences, Harvard University, United States


In biological contexts as diverse as development, apoptosis, and synthetic microbial consortia, collections of cells or subcellular components have been shown to overcome the slow signaling speed of simple diffusion by utilizing diffusive relays, in which the presence of one type of diffusible signaling molecule triggers participation in the emission of the same type of molecule. This collective effect gives rise to fast-traveling diffusive waves. Here, in the context of cell signaling, we show that system dimensionality – the shape of the extracellular medium and the distribution of cells within it – can dramatically affect the wave dynamics, but that these dynamics are insensitive to details of cellular activation. As an example, we show that neutrophil swarming experiments exhibit dynamical signatures consistent with the proposed signaling motif. We further show that cell signaling relays generate much steeper concentration profiles than does simple diffusion, which may facilitate neutrophil chemotaxis.


Prototypical diffusive signaling – in which individual cells communicate with neighbors by releasing diffusible molecules into the extracellular medium – is a relatively slow process. Signaling molecules undergoing random walks in the extracellular medium have a root mean square displacement that grows like the square root of both the time since emission, t, and the signaling molecule diffusivity, D. It follows that the distance an individual cell can signal also grows like the square root of time. Thus, for thousands of cells coordinating actions over millimeters, simple diffusive signaling with small molecules (D10-10 m2/s) takes hours. These length and times scales are incommensurate with observed behavior in developmental biology (Chang and Ferrell, 2013; Cheng and Ferrell, 2018; Vergassola et al., 2018), immune response (Reátegui et al., 2017), and microbial consortia (Parkin and Murray, 2018), in which cells exchanging diffusible molecules coordinate activity over millimeters in tens of minutes.

Indeed, when many cells collectively integrate environmental cues and participate in the signaling, they can propagate diffusive waves with a fixed speed, v, in the asymptotic limit. This effect and its analogs have long been studied in the context of excitable media (Keener and Sneyd, 2009; Keener, 1987; Muratov, 2000) and observed in biological phenomena as diverse as natural cell signaling circuits (Noorbakhsh et al., 2015; Pálsson and Cox, 1996; Kessler and Levine, 1993; Gelens et al., 2014), synthetic cell signaling circuits (Parkin and Murray, 2018), apoptosis (Cheng and Ferrell, 2018), range expansions (Tanaka et al., 2017; Fisher, 1937; Kolmogorov et al., 1937; Barton and Turelli, 2011; Gandhi et al., 2016; Birzu et al., 2018), and development (Chang and Ferrell, 2013; Vergassola et al., 2018; Muratov and Shvartsman, 2004; Nolet et al., 2020). In this way, small groups of cells can transmit signals more quickly than simple diffusion allows by recruiting the help of their neighbors.

While diffusive waves have been observed in a variety of biological processes, they have also been experimentally probed in a variety of spatial contexts – in quasi-1D tubes (Cheng and Ferrell, 2018; Chang and Ferrell, 2013; Nolet et al., 2020), in quasi-2D droplets and chambers (Nolet et al., 2020; Afanzar et al., 2020), on 2D surfaces in fly eggs (Vergassola et al., 2018), and on substrates of finite thickness (Parkin and Murray, 2018; Pálsson and Cox, 1996). And while the phenomenology of diffusive waves has been studied for years, in the context of cell signaling it is less well-understood how the propagation and initiation of such waves are affected by the dimensionalities of the cellular distribution and the diffusive environment – or even how to identify the system dimensionality – as previous modeling work has largely assumed quasi-1D dynamics (Kessler and Levine, 1993; Meyer, 1991; Gelens et al., 2014; Vergassola et al., 2018). Also unclear is how robust the resulting signaling dynamics are to underlying biological details, such as the shape of the function governing cell activation and signaling molecule emission.

Here, we revisit the propagation and initiation of diffusive waves in the context of cell signaling. Through a comprehensive study of single-component relays — in which cells measure the local concentration of a signaling molecule and participate in the emission of the same molecule — we show that the asymptotic wave dynamics of diffusive relays are governed by simple scaling laws. In some system dimensionalities, these scaling laws are identical to famous results from the 20th century (Fisher, 1937; Kolmogorov et al., 1937; Luther, 1906); in other system dimensionalities, we show that these well-known scaling laws can be drastically altered. For example, cells confined to two (or one) dimensions with signaling molecule diffusion in three (or two) dimensions give rise to a diffusive wave whose speed has no dependence on D: a wave driven by diffusion whose speed does not depend on the rate of diffusion. In contrast to the dramatic effect of system dimensionality, these scaling laws are insensitive to many biological details, including the functional form of cellular activation — the dependence of signaling molecule emission rate on the local concentration. We additionally account for other phenomena – molecule decay, pulsed emission, and the discreteness of cells – that do affect the asymptotic wave dynamics; in so doing, we provide an intuitive rubric for determining under what conditions these effects alter the wave propagation speed.

In our studies of wave initiation, we systematically examine under what conditions a group of cells can trigger the formation of a diffusive wave. Here again, our results provide predictive relationships between biophysical inputs and the resulting dynamics, which are at once dramatically affected by dimensionality and largely insensitive to the details of activation and cellular uptake.

Finally, we show that neutrophil swarming experiments (Reátegui et al., 2017) display dynamics consistent with our model. In this context, our results elucidate a potential design principle of diffusive relays: they create large concentration gradients. Whereas simple diffusion of a signaling molecule from a central source creates a shallow concentration profile that falls off like exp(-r2/4Dt), relays give rise to steep concentration profiles with gradients that quickly propagate outward and decay only modestly inside the wave front. As such, for cells like neutrophils – which use a small molecule, leukotriene B4 (LTB4), as an intercellular signaling molecule and chemoattractant (4, 18, 19) – relays may provide a method for cells to collectively generate large, continous chemical gradients that may serve to guide directional migration; the continuous gradients generated by single-component relays contrast with the pulse trains of chemotactic cues observed in, for example, Dictyostelium discoideum (Kessler and Levine, 1993; Pálsson and Cox, 1996).


Model construction

We begin by considering a static group of cells uniformly distributed in two dimensions – for example, atop a solid surface – and described by an area density ρ (Figure 1A). We assume a cell at position 𝐫 senses the local concentration of a signaling molecule, c(𝐫,t), and participates in the emission at a concentration-dependent rate af(c) with a the maximum rate and f(c) a dimensionless function. Once secreted into the extracellular medium, the signaling molecules diffuse with diffusivity D. Treating the cells and signaling molecule concentration in the continuum limit – we discuss the validity of doing in the next section and in Appendix 6: Assessing the validity of a continuum analysis – gives rise to a single equation that governs the time evolution of c(𝐫,t):

(1) ct=D2c+aρδ(z)f(c)

where the Dirac delta function δ(z) accounts for the fact that the cells are confined to the plane. The source function f(c) is in general a complicated non-linear function of c. It can include uptake, release, and cell-induced degradation of the signaling molecule – or any other process proportional to the local cell density. We will consider this general case shortly. To start, we consider a simple case in which cells measure the local signaling molecule concentration, c, and participate in the emission only if c exceeds a threshold concentration, Cth. In such a case, the activation function f(c) is well described by a Heaviside step function Θ[c-Cth] and the concentration dynamics obey

(2) ct=D2c+aρδ(z)Θ[c-Cth].

Additionally, while we at first consider cells scattered in a two-dimensional plane, one can study the signaling dynamics of cells in a one-dimensional channel or a three-dimensional environment with similar analyses. Below, we discuss the connections between the cell signaling dynamics in all these scenarios, and all are treated in depth in Appendix 2: Asymptotic wave ansatz.

Asymptotic relay dynamics of cells in 2D with diffusion in 2D.

(A) Schematic illustrating the diffusive relay motif. Cells (pink with purple nucleus) release a signaling molecule that diffuses (blue clouds). They do so when the local concentration exceeds a threshold, Cth. This gives rise to a diffusive wave with wave speed v. (B) Snapshot concentration profiles. Asymptotic theory (Equation (6), black lines) and numerical simulation of Equation (2) (red dots, details of the numerical methods can be found in Materials and methods) are in good agreement and show outward-propagating waves. Here, D=10-10 m2/s, v=2 µm/s, and hCth/aρ=D/v2. Numerical simulations assume that a cell colony of size ri=4D/v (dashed vertical line) centered at the origin starts signaling t=0. (C) Numerical wave speed as measured at t=100D/v2 (markers) agrees well with theory (Equation (5), black line) as we independently vary D (circles) and hCth/aρ (diamonds) relative to the panel B values (red circle and diamond).

Asymptotic wave dynamics

Our first step in understanding diffusive signaling relays is to solve for the asymptotic dynamics of Equation (2). Since such relays involve cells signaling their neighbors, which then signal their own neighbors, one can imagine that diffusive relays give rise to diffusive waves. We therefore make the ansatz that the concentration c(𝐫,t)=c(r,z,t) can be described by an outward-traveling wave of the form c(r,z,t)=c(r~=r-vt,z) (Fisher, 1937; Kolmogorov et al., 1937; Tanaka et al., 2017). Here, r~ is the distance from the wave front – negative when inside the wave front, positive when beyond – and v is the wave speed. In essence, we wish to examine the wave from the perspective of an observer moving at the wave front. With Cthc(r~=0,z=0) and rD/v, we take Equation (2) and arrive at the following equation governing asymptotic behavior:

(3) 0=D(2cr~2+1rcr~+2cz2)+vcr~+aρδ(z)Θ[cCth]D(2cr~2+2cz2)+vcr~+aρδ(z)Θ[cCth]=D(2cr~2+2cz2)+vcr~+aρδ(z)Θ[r~].

Since we consider rD/v, we may ignore the D(c/r~)/r term due to the dominance of vc/r~. This is effectively the same as ignoring the curvature of the wave front and has the effect of reducing our asymptotic analysis of cells in two dimensions into an asymptotic analysis of cells in one dimension (Tanaka et al., 2017). The asymptotic dynamics of cells distributed in three spatial dimensions allow for a similar manipulation (see Appendix 2: Asymptotic wave ansatz).

We wish to find a solution to Equation (3) for various diffusive – that is, extracellular – environments. In doing so, we hope to solve for the spatial dependence of the concentration profiles c(r~,z) as well as a relationship that will tell us how the signaling dynamics – in this case, the wave speed v – depend on the biophysical system parameters like the cell density, ρ; the concentration threshold, Cth; and the signaling molecule emission rate, a.

But first, we note that Equation (3) provides two quantities of value: a natural length scale D/v and a natural time scale D/v2. For a small diffusing molecule with D10-10 m2/s and a wave speed of v1 µm/s – approximately the numbers relevant for several experimental systems (Cheng and Ferrell, 2018; Chang and Ferrell, 2013; Parkin and Murray, 2018; Vergassola et al., 2018; Pálsson and Cox, 1996) including, as we show below, neutrophil swarming (Reátegui et al., 2017) – we recover D/v100 µm and D/v2100 s. We have already used the natural length scale D/v to derive Equation (3) and to show that cells in 2D have the same asymptotic dynamics as cells in 1D or 3D, and we can use these scales to further justify several other approximations we have made so far. For instance, the approximation that the out-of-plane cell density can be described by δ(z) is valid when the cell size HD/v; similarly, decay of the signaling molecule can be neglected for a decay rate γ(D/v2)-1 while pulsed emission gives rise to the same asymptotic wave speed if the width of the pulse τ satisfies τD/v2. Finally, we note that the use of Equation (2) as a starting point is justified when the mean distance d between neighboring cells satisfies dv/4D1. A thorough, mathematical discussion of all the above, including a demonstration of why D/v and D/v2 are the appropriate scales, is presented in Appendix 2: Asymptotic wave ansatz.

When the extracellular medium thickness hD/v, diffusion of the signaling molecule is effectively two-dimensional as we can take 2c/z20 and δ(z)1/h. In this limit, Equation (3) becomes

(4) hD/v:0=D2cr~2+vcr~+aρhΘ[cCth]=D2cr~2+vcr~+aρhΘ[r~]

which we can solve to find the asymptotic dynamics of cells in 2D (1D, 3D) with effective signaling molecule diffusion in 2D (1D, 3D) – the thin extracellular medium limit (Figure 1). This corresponds to the long-pulse, long-decay time limit of the model constructed by Kessler and Levine, 1993 and is similar to the model considered by Meyer, 1991. Adding signaling molecule decay to Equation (4) would yield a model first considered by McKean, 1970 in the context of nerve impulse propagation.

Before solving Equation (4) exactly, we make two crucial observations from which we can derive the functional form of the wave speed, v. First, because the source (furthest right) term in Equation (4) is proportional to aρ/h, all concentrations in the problem, including Cth, are proportional to aρ/h. As the non-source terms in Equations (4) and (1) are linear, the only role aρ/h serves is to set the concentration scale of the dynamics. Thus, Cth, a, ρ, and h combine to give us a single model parameter to describe the threshold concentration, hCth/aρ, which has units of time (measured in s). Second, the only other parameter in the problem besides v – which we want to calculate – is the diffusion constant, D, which has units of length squared divided by time (measured in m2/s). Thus, the only combination of these two parameters that will give a speed (measured in m/s) is (aρD/hCth)1/2. By this simple dimensional analysis argument, the wave speed v can only be v=α(aρD/hCth)1/2 for some constant α. Formally, the above procedure is equivalent to non-dimensionalizing Equation (4), as discussed in Appendix 2: Asymptotic wave ansatz.

By the same reasoning, any activation function f(c) – a Heaviside step function, a Hill function, or even a bistable function – that can parameterized by a single concentration Cth and emission rate a must give the same scalings if it has a traveling wave solution. While we focus on positive activation functions in this work, we emphasize that if signaling molecule degradation is dominated by cell-induced processes like uptake, then signaling molecule degradation is also proportional to the cell density and the resulting (presumably bistable) production curve will yield dynamics that are also beholden to this scaling law.

One can confirm this scaling law for Heaviside activation by solving Equation (4) for r~>0 and r~<0, then matching boundary conditions at r~=0. This analysis indeed reveals that

(5) hD/v: Cth=aρD/hv2  v=aρD/hCth


(6) hD/v:c(r~)={aρr~/hv+aρD/hv2r~0aρDer~v/D/hv2r~0.

The concentration of signaling molecule thus grows linearly in the distance inside the wave front and decays exponentially in the distance beyond the wave front. We compare numerical simulations of Equation (2) (see Materials and methods for details) with the above asymptotic formulae for wave speeds and concentration profiles in Figure 1B/C. For rD/v, the asymptotic formulae describe well both the concentration profile and the wave speed.

The wave speed relationship given in Equation (5) is analogous to the Fisher-Kolmogorov wave speed (Fisher, 1937; Kolmogorov et al., 1937; Gelens et al., 2014) – with hCth/aρ replacing the doubling time as the characteristic time scale in the problem – and has been discussed in beautiful previous work (Gelens et al., 2014; Meyer, 1991), starting with Luther, 1906. Amazingly, Luther’s formula, which posits the scaling relation vD, holds even in scenarios beyond those considered here; for instance, waves driven by oscillatory activation dynamics – as are relevant for intercellular signaling in Dictyostelium discoideum (Kessler and Levine, 1993; Pálsson and Cox, 1996) and developmental trigger wave propagation (Gelens et al., 2014; Chang and Ferrell, 2013) – are subject to this same scaling. One can understand this through simple dimensional analysis. These more complex scenarios add signaling molecule decay and a periodically modulated source function to the above model. Thus, to our set of parameters, D (measured in m2/s) and hCth/aρ (measured in s), we add a modulation time τ (measured in seconds) and decay rate γ (measured in 1/s). As D is the only parameter involving a length scale, it must be that vD even in these more complex scenarios.

By way of contrast, Vergassola et al., 2018 have shown that an unconventional scaling of vD3/4 can result from time-dependent dynamics of the source term at the wave front, a phenomenon that breaks our assumption that all cells obey the same time-independent source function f(c). Similarly, as we will now show, the dimensionality of the system can also have a dramatic effect on wave speed scaling laws.

Next, we consider a thick extracellular medium for which hD/v. Such a configuration is relevant for signaling in bacterial consortia atop thick, permeable substrates (Parkin and Murray, 2018) or anywhere that a lower dimensional tissue abuts a thick and permeable extracellular environment as can be found, for example, in the retina. Here, the signaling molecules can diffuse out of plane (Figure 2A). Because the cells sit atop a solid boundary, signaling molecules can only diffuse in the upper half of the plane and the source term in Equation (3) acquires a factor of two to account for this boundary condition:

(7) hD/v:0=D(2cr~2+2cz2)+vcr~+2aρδ(z)Θ[cCth]=D(2cr~2+2cz2)+vcr~+2aρδ(z)Θ[r~].
Asymptotic relay dynamics with cells in 2D and diffusion in 3D.

(A) Schematic of cells (pink with purple nucleus) performing a diffusive relay in which signaling molecules (blue clouds) can diffuse out-of-plane. Here, such relays give rise to a diffusion-constant-independent wave speed, v. (B) Snapshot concentration profiles of the signaling molecule show good agreement between numerical simulation of Equation (2) (blue dots, details of numerical methods can be found in Materials and methods) and asymptotic theory (Equation (9), black lines). Here, D=10-10 m2/s and v=2 µm/s with Cth/aρ=2/πv. The initial signaling colony is of size ri=4D/v (dashed vertical line). (C) Numerical wave speed as measured at t=100D/v2 (markers) agrees well with theory (Equation (8), black line) as we independently vary D (circles) and Cth/aρ (diamonds) relative to the panel B values (blue circle and diamond). As predicted, v is indeed D-independent in this system.

Effectively, we have cells in 2D with diffusion in 3D. We note that this case is asymptotically equivalent to cells in 1D emitting into a semi-infinite 2D environment. Thus, comparing to Equation (6), we can see that the asymptotic dynamics are not determined by the dimension of the cell distribution or the diffusive environment, but by the difference in dimension between them.

The same dynamics hold for cells on a curved surface (such as epithelia) as long as the length scale of the curvature and the thickness of the extracellular medium are both large compared to D/v. If the length scale of the curvature is large compared to D/v, but the extracellular medium is thin compared to D/v, then the dynamics will be of cells in a 2D plane with diffusion in 2D. Similarly, cells on the surface of a tube with diffusion in the tube's interior will interpolate between these two limits: when the radius of the tube is large compared to D/v, the dynamics will be of cells in 2D and diffusion in 3D; when the radius of the tube is small compared to D/v, the dynamics will be of diffusion and cells in 1D.

Examining Equation (7) as we did Equation (4) reveals that every concentration in a thick extracellular medium is proportional to aρ. Thus, we have two independent parameters in Equation (7): Cth/aρ (measured in s/m) and D (m2/s). The only combination of these parameters that will give a wave speed (measured in m/s) is aρ/Cth. It therefore must be the case that v=αaρ/Cth with α a constant – a wave driven by diffusion whose wave speed is independent of the rate of diffusion. We again stress that this is true for any activation function that has a traveling wave solution and can be parameterized by a single concentration Cth and a single emission rate a. (For Hill function activation, α2/π for n2, see Appendix 5: Asymptotic wave dynamics with Hill function activation.) Thus, the scaling laws governing the asymptotic dynamics are insensitive to the details of single-cell activation.

It is worth reflecting on the fact that some system geometries give a wave whose speed is diffusion constant-independent. This finding implies that, at least in some contexts, the size of the signaling molecule has little to do with the resultant cell signaling speed. We note that this is in contrast with the more standard wave speed scaling in Equation (5), in which smaller (lower molecular weight, higher D) signaling molecules result in a faster wave, all else equal.

A full solution of Equation (7), obtained in Appendix 2: Asymptotic wave ansatz by combining a partial Fourier transform in the z-dimension and the methods used to solve Equation (4), yields

(8) hD/v:Cth=2aρ/πvv=2aρ/πCth


(9) hD/v: c(r~){2aρ(r~/πvD)1/2r~D/vaρ(D/πr~v3)1/2evr~/Dr~D/v.

So for cells in a thick extracellular medium, the concentration grows like the square root of the distance inside the wave front and decays exponentially beyond the wave front. As with the 2D diffusive environment, we verify these relationships numerically (Figure 2B/C, see Materials and methods for details of the numerical simulations). We see that the wave speed is indeed D-independent over two orders of magnitude in the diffusion constant.

The diffusive relay signaling motif therefore gives rise to diffusive information waves for which Equation (5) and Equation (8) provide predictive relationships between wave speed, threshold concentration, cell density, extracellular medium thickness, and emission rate for a variety of system dimensionalities. Similarly, Equation (6) and Equation (9) provide quantitative functional predictions of the concentration profiles generated by diffusive relays. By dimensional analysis, these scaling laws are insensitive to the details of activation. Nonetheless, other details – signaling molecule decay, pulsed emission, discreteness of cells – can alter these robust scaling laws (Keener, 2000; Dieterle P and Amir A, 2020. Manuscript in preparation). We explicitly discuss these corrections in the appendices (see Appendix 3: Pulsed emission and decay and Appendix 6: Assessing the validity of a continuum analysis), where we also discuss the dynamics of cells in 1D with 3D diffusion and the properties of waves in an arbitrary extracellular medium thickness. Both signaling molecule decay and pulsed emission decrease the steepness of the concentration gradient inside the wave front, and both decrease the wave speed. We emphasize that, in all cases, the asymptotic dynamics are not determined by the dimension of the diffusive or cellular environment, but by the difference in dimension between the two.

Signaling wave initiation

Armed with a knowledge that diffusive relays birth diffusive waves, we now ask whether such waves are always initiated. As with the asymptotic dynamics, wave initiation depends on the system dimensionality. Here, however, the dimensionality of the diffusive environment alone determines qualitative behavior. Much previous work in chemical waves and excitable media has shown that a delicate interplay of activation, repression, and diffusion can give rise to a host of dimension-dependent wave initiation phenomena (Foerster et al., 1989; Weise and Panfilov, 2011); our task here is to study the dimension-dependent dynamics of concentration build up in single-component relays.

To begin, we consider an ‘initiating colony’ of radius ri in which cells emit a diffusible signaling molecule with rate a (Figure 3A). The surrounding cells respond by emitting the same signaling molecule according to some activation function, f(c). Here, we take f(c) to be a Hill function of degree n (Figure 3B).

Wave initiation dynamics.

(A) Schematic demonstrating wave initiation. Cells within some initial signaling volume of radius ri begin signaling at some rate a. The signaling wave is initiated when the concentration at nearby cells exceeds the threshold concentration, Cth. (B) Cells near the initial signaling volume participate in the emission according an activation function, f(c). For instance, in the case of Hill function activation, f(c)=[1+(Cth/c)n]-1. C: Initiation times for Heaviside activation, in which f(c)=Θ[c-Cth]. Numerics (thick colored lines) and approximate asymptotic theory (Equations (10) and (11), thin black lines) of the initiation time’s dependence on ri for cells and diffusion in 1D or 2D (left) or cells in 2D with diffusion in 3D (right). For cells and diffusion in 1D, Equation (10) provides a good approximation in the limits vri/D1 and vri/D1. Similarly, for cells and diffusion in 2D, Equation (11) governs the large and small vri/D limits. In both of these cases, the wave always initiates, but the initiation time can be orders of magnitude larger than D/v2 if riD/v. For cells and diffusion in 3D (right), signaling waves do not initiate for vri/D<3. Here again, the asymptotic theory Equation (12) is in good agreement with numerics.

In one- and two-dimensional diffusive environments, a continuously emitting source leads to a diverging concentration throughout the space. However, in three-dimensional diffusive environments, a continuously emitting source gives rise to a steady-state concentration with 1/r tails (Krapivsky et al., 2010).

We therefore expect that an initiating colony of cells, regardless of its radius ri (in fact, even if it consists of a single cell – see Appendix 7: Initiation dynamics), will be able initiate a diffusive wave in one- and two-dimensional environments; meanwhile, a colony in a three-dimensional environment may fail to initiate a diffusive wave. This is indeed what we observe.

In Figure 3C, we demonstrate this dramatic dimension-dependence in the case of switch-like activation, for which f(c)=Θ[c-Cth]. To find the initiation time tinit, we integrate Green’s functions of the diffusion equation (see Appendix 7: Initiation dynamics for details) to calculate the concentration profile created by cells continuously emitting with rate a inside the initiating colony. When the concentration at ri is equal to the threshold – Cth=c(ri,tinit) – cells outside the initiating colony begin to participate in the relay and the wave is initiated. Below, we characterize the initiation time as a function of ri and the characteristic time and length scales – D/v2 and D/v, respectively – of a given system, thus linking the initiation dynamics to the asymptotic wave speed, v. A summary of these results as a function of dimension, along with a summary of the asymptotic dynamics, can be found in Table 1.

Table 1
Summary of asymptotic and initiation dynamics with Heaviside activation.

For different system dimensionalities, we summarize the asymptotic wave speed, v; the initiation time for small initial signaling colony size, tinit,vriD1; and the initiation time for large initial signaling colony size, tinit,vriD1. One-dimensional diffusive environments are assumed to be narrow channels of width h in each direction perpendicular to the channel length. The cell density ρ has units 1/m for cells in 1D, 1/m2 for cells in 2D, and 1/m3 for cells in 3D. When the diffusive and cell dimensions do not match, the environment is assumed to be semi-infinite.

Cells in 1D, diff. in 1D(aρDh2Cth)1/2(Dvri)22D/v2
Cells in 1D, diff. in 2D2aρπhCthexp(2Dvri)24D/πv2
Cells in 2D, diff. in 2D(aρDhCth)1/2exp(2Dvri)22D/v2
Cells in 2D, diff. in 3D2aρπvno waves4D/πv2
Cells in 3D, diff. in 3D(aρDCth)1/2no waves2D/v2

For cells in 1D with 1D diffusion, the initiation time in the limits of small (riD/v) and large (riD/v) initiating colonies is:

(10) tinit{(πD/4v2)(D/vri)2riD/vtmin,1,1=2D/v2riD/v.

When riD/v, the signaling molecules quickly diffuse across and away from the initiating colony. Thus, it is hard for the colony to build up a concentration that exceeds Cth. Correspondingly, the initiation time increases like 1/ri2 for small ri. Meanwhile, for riD/v, the size of the initiating colony becomes irrelevant and reaches a minimum value of tmin,1D, determined entirely by the characteristic time scale D/v2. The full dependence of tinit on ri is pictured in Figure 3C, where we show that the above limits are valid approximations.

Next, we consider cells in 2D with diffusion in 2D. Here, for riD/v, the initiation time scales harshly as

(11) tinit{(ri2/4D)e(2D/vri)2riD/vtmin,2,2=tmin,1,1=2D/v2riD/v.

Results in the above limits are corroborated by numerical simulation in Figure 3C, where we show initiation times for the limits above and for intermediate values of vri/D.

Lastly, we consider cells in 3D with 3D diffusion and find that there is a critical initial signaling colony size of ri=3D/v below which the wave will not initiate. Around the critical colony size, tinit diverges as (vri/D)6[(vri/3D)2-1]-2. If riD/v, then tinit again plateaus at a constant value tmin,3D that only depends on the characteristic time scale D/v2:

(12) tinit{no initiationri<3 D/v(D/9πv2)(vri/D)6[(vri/3D)21]2ri3 D/vtmin,3D=2D/v2ri3 D/v.

These analytic expressions again agree well with numerical simulation, as seen in Figure 3C. In Appendix 7: Initiation dynamics, we work out the case of cells in 2D with diffusion in 3D, for which there is a minimum initiating colony size of ri=D/v. There, we also show that the qualitative findings presented above also hold for systems with discrete cells.

The critical initiating colony size for a 3D environment is reminiscent of elegant work on range expansions (Tanaka et al., 2017; Barton and Turelli, 2011). There, the effects of diffusive migration and population growth compete with each other, and a critical mass is needed to initiate the spatial advance of a particular genotype. Here, the dimension-dependent dynamics of concentration build-up dictate that a signaling wave which will always initiate in one- and two-dimensional environments requires a critical initial colony size in 3D.

Because the signaling wave always initiates in one- and two-dimensional environments, it can in principle be initiated by a single cell. As random activation of a single cell can initiate a signaling wave that fixes the entire population to maximal activation, these signaling dynamics have typically been thought of as unstable (Deneke and Di Talia, 2018). Yet, as we have shown here, even in one- and two-dimensional environments, the initiation time for colonies smaller than D/v can be many orders of magnitude larger than the characteristic time scale of D/v2 (Figure 3D). Thus, even though this signaling modality is technically unstable, it is robust against stochastic activation of a small number of cells over very long time scales.

In effect, then, even strictly positive-valued activation functions require a ‘critical mass’ of cells to initiate a signaling wave. In the context of neutrophil swarming – which we will shortly consider in more detail – this critical mass may provide a basis by which the immune system ‘decides’ whether to initiate a full-scale swarming response. In vitro experiments (Reátegui et al., 2017) indicate that small colonies of a pathogen can indeed fail to incite a swarm. Moreover, since the critical size of an initiating colony goes like D/v, we can see that relays utilizing smaller (lower molecular weight, higher D) signaling molecules require larger critical masses, all else equal.

Finally, we note that for cells with a Hill-like activation function f(c)=cn/(cn+Cthn) of order n2, the above results for switch-like activation provide a good quantitative approximation of the initiation times (see Appendix 8: Wave initiation with Hill function activation). Moreover, for cells in 3D with Hill activation functions of order n>3, there is a critical colony size just as for switch-like activation. These results highlight the role of spatial degrees of freedom in determining the wave initiation dynamics and stability.

Application to neutrophil swarming and gradient generation

With a firm understanding of the diffusive wave and initiation dynamics, we now turn our sights to understanding a specific model system: neutrophil swarming. In beautiful work across several organisms (Lämmermann et al., 2013; Isles et al., 2019; Reátegui et al., 2017), experimentalists have observed striking behavior: an acute injury or infection can elicit rapid, highly directive motion of neutrophils – the most prevalent white blood cells – toward the site of the injury or infection. These experiments have demonstrated that a lipid small molecule called leukotriene B4 (LTB4) – along with many larger, slower-diffusing proteins (Reátegui et al., 2017) – governs the long-range recruitment of swarming neutrophils (Lämmermann et al., 2013; Afonso et al., 2012; Reátegui et al., 2017; Isles et al., 2019). Reategui et al. have noted the presence of several other pro- and anti-inflammatory lipid small molecules during swarming, though their precise roles are less clear. LTB4 serves to activate the neutrophils and also acts as a chemoattractant (Afonso et al., 2012) when receptors for LTB4 are blocked, swarming behavior is significantly impaired (Lämmermann et al., 2013; Reátegui et al., 2017). The release of LTB4 has been thought to work as a relay, although the precise mechanistic details of this relay remain unclear (Lämmermann and Germain, 2014; Kienle and Lämmermann, 2016).

In vitro experiments performed with human neutrophils are particularly relevant given the results discussed so far. In these experiments, human-derived neutrophils are injected into a chamber, then settle onto the surface of a glass slide, resulting in a uniform sprinkling of cells in 2D. Also on the glass slide are circular 'targets' (of size ri) coated in zymosan, a fungal surface protein that elicits a swarming response (Reátegui et al., 2017). Some cells land on or near the target, giving an initial condition as in Figure 3A. These cells begin signaling their neighbors, which in turn migrate towards the target (Figure 4A).

Application to neutrophil swarming.

(A) Schematic of the simple diffusion model. Here, cells on the target (within ri) signal distant neighbors by continuously emitting a single signaling molecule. If the neighboring cells have a chemotactic response, they migrate toward the target with some noise – that is, some non-zero angle θ with respect to the target. Otherwise, they move around with no sustained directionality. (B) Experimental data (color plot) reproduced from Reátegui et al., 2017 showing the information wave front in neutrophil swarming experiments. By tracking the neutrophils in space and time, they observe highly directed motion of the neutrophils towards the target (pink) starting around t=200 s. There is a clear boundary in space and time – the information wave front – between the regions where cells migrate toward the target (pink) and jostle around with no particular direction (white and light blue). While a relay theory (black line) is consistent with the convex shape of the information wave front, simple diffusive signaling by only the cells on the target (gray line) is not. The diffusion constants for both models is D=1.25×10-10 m2/s. The threshold concentrations for the relay and simple diffusion models are Cth/aρ3.66×105 s/m and 2.91×104 s/m, respectively. The parameters for the relay model are chosen to fit the wave front by eye while the simple diffusion model parameters are chosen to give the same signaling distance at t=500 s. (C) Gradients created by signaling relays (black) and simple diffusion (gray) models in panel B. The dashed vertical lines indicate the location of the information wave front. As time increases from left to right, the relay signaling motif gives an information wave that signals cells faster than simple diffusion in the long time limit. Cells within the wave front (to the left of the dashed lines that indicate the wave fronts) experience significantly larger gradients when the cells utilize a relay, which may facilitate efficient chemotaxis.

By tracking individual cells in time, one can deduce their migratory direction as a function of time. A typical metric for quantifying the directionality a cell’s migration is the chemotactic index – the cosine of the angle θ between a cell’s motion and the direction of the target (Figure 4A). One can average over the cells at a given distance r and time t to construct a plot of the average directionality cosθ in space and time. As pictured in Figure 4B, such a plot reveals a clear divide in space and time between cells that are highly directed toward the target (pink) and those without any particular directionality (white and light blue). We refer to the boundary of this divide as an information wave front – cells that lie underneath the curve have received the signal and begun chemotaxing toward the target while those above the curve have not.

Interestingly, the information wave front is convex with respect to the origin – a dramatic departure from what simple diffusive signaling by cells on the target would yield (Figure 4A/B), and from what Reategui et al. observe in experiments with neutrophils whose LTB4 receptors have been blocked (see Appendix 10: Simple diffusion model for more). We therefore posit that the cells may be participating in a relay in which they emit LTB4 in response to the same and check to see if this is consistent with the observed information wave front.

To do so, we perform a numerical simulation of Equation (2) with an additional term to account for the signaling of cells that land on the target. For this analysis, we assume a circular target of radius ri100 µm, though the targets fabricated by Reategui et al. are smaller, oblong objects. Here, the diffusive environment is effectively three dimensional and the cells are close enough to allow for the use of a continuum model like Equation (2) (see below). Our model assumes switch-like activation of neutrophils, which we associate with the onset of directed chemotaxis. We ignore the inward migration of cells in this analysis, as it has a negligible effect on the information wave propagation since the cells move at a speed u0.3 µm/s v (see Appendix 11: Quantifying the effects of chemotaxis). Thus, as mentioned above, Equation (2) effectively has two parameters: Cth/aρ and D. Fitting these two parameters to the observed information wave front gives Cth/aρ3.67×105 s/m and D1.25×10-10 m2/s, the latter of which is consistent with the diffusion constant of a small molecule like LTB4. This implies a wave speed of v1.7 µm/s. Thus, we are validated in using a continuum model with a thick extracellular medium, as for this experiment the extracellular medium thickness h=2 mm D/v and the mean distance between neutrophils, d=50 µm, satisfies vd/4D0.171. The cell thickness H10 µm indeed satisfies HD/v, meaning the use of the delta function to describe the cell distribution is valid. Finally, as LTB4 has a lifetime 1/γ of many minutes (Bray, 1983) and D/v240 s 1/γ, we can indeed ignore signaling molecule decay. These fit parameters give a curve that matches the transient dynamics over the field of view of the experiment (Figure 4). Thus, our relay model gives dynamics that are consistent with the dynamics of neutrophil swarming experiments – namely, the observed convex shape of the information wave front. Larger field-of-view and longer time-course experiments with varying cell densities and larger targets will provide a deeper mechanistic understanding of such relays, while also testing the scaling predictions of Equation (5) and Equation (8).

The fit value of Cth/aρ=3.67×105 s/m is consistent with the neutrophil’s LTB4 receptor affinity. To show this, we first note that Reategui et al. measured the LTB4 emission rate under similar conditions as the relay experiment analyzed above; they found that a40 molecules per second per cell (see Appendix 9: Sensitivity of the information front to fit parameters for details). Using the cell density of ρ=1/d2=(50μm)-2, we find that Cth500 pM. This value is within the range of the measured BLT1 receptor affinity for LTB4, which is reported to be approximately 0.1 − 2 nM (Yokomizo, 2015).

Finally, we comment on the matter of why neutrophils might employ such signaling relays. As we have shown above, relays lead to 'fast' communication, in the sense that they give rise to diffusive waves which travel a distance vt in a time t, compared to the Dt distance of simple diffusion. However, there is another potential reason to use diffusive relays: they create strong gradients that may help cells chemotax effectively.

To get an idea of the gradients we are working with, we compare those generated by a relay – calculated by solving Equation (2) and approximated in Equation (9) – to a comparable simple diffusion model, such as that pictured in Figure 4B. (In Appendix 10: Simple diffusion model, we present the same comparison for a thin extracellular medium.) As is well-known, a burst-like emission of a diffusible molecule creates shallow, Gaussian concentration profiles away from the source; the same is true for continuous emission of a fixed source. Thus, the gradients that individual cells or small colonies of cells can create through simple diffusive signaling are orders of magnitude shallower than the collective gradients generated by relays (Figure 4C). This hints that neutrophils may use relays not solely for their improved signaling speed, but also for the strong resulting chemotactic gradients.


In this work, we have shown how simple cell signaling relays can give rise to diffusive waves whose properties are robust to many underlying details. Our work especially highlights the importance of the dimensionality of the extracellular medium, as seemingly innocent changes to the environment can have large effects on the resulting diffusive waves. The strong effect of system dimensionality is reminiscent of previous work on diffusive dynamics, which showed how dimensionality can effect Turing pattern instabilities (Levine and Rappel, 2005).

Although we have characterized the asymptotic dynamics, initiation, and potential design principles of these waves in several scenarios, many interesting problems remain as yet unsolved. First, as noted by Lammermann and colleagues (Lämmermann and Germain, 2014; Kienle and Lämmermann, 2016), it is unclear how the complexities of in vivo extracellular environments affect these results, particularly in the context of neutrophil swarming. Ambient flow (for example, in blood vessels), constrictions, and complex diffusive environments may lead to dynamics of biological relevance beyond those discussed here. Additionally, it would be interesting to study how different models of chemoreception and cellular uptake – topics of theoretical (Muratov and Shvartsman, 2004) and experimental (Youk and Lim, 2014; Scherber et al., 2012; Tweedy et al., 2016) relevance – affect our conclusions.

As an experimental test of our model, we propose studying neutrophil swarming dynamics over a wide field of view with varying cell densities and extracellular medium thicknesses. For diffusive waves with approximately our experimentally inferred parameters for neutrophil swarming (D/v100μ µm), one could probe the thin extracellular medium limit of hD/v with microfluidic chambers of tens of microns in height. Similarly, with mm-scale chambers introduced by Reátegui et al., 2017 and discussed in the previous section, one can reach the limit of a thick extracellular medium. Experiments in these two limits would provide quantitative tests of our theory. In particular, varying cell density would provide a test of the dimensionality-dependent relations for collective signaling wave speed, Equations (5) and (8).

On a mechanistic level, although a relay mechanism would allow neutrophils to quickly coordinate their response, it remains unclear how inflammatory response is modulated in such a scenario. If inflammation during neutrophil swarming is governed by a fast-travelling wave, then how do the cells collectively turn off response? One possibility is that signaling pathways in neutrophil swarm resolution – for instance, those involving LXA4 (Reátegui et al., 2017) production and emission – work by a similar relay mechanism; it is also possible that LTB4 production is governed by other fast-diffusing signaling molecules whose presence is necessary for LTB4 production, thereby limiting the relay’s recruitment range.

Studies of the neutrophil relay mechanism may provide an interesting contrast to similar intercellular signaling dynamics in Dictyostelium discoideum (Pálsson and Cox, 1996; Kessler and Levine, 1993; Noorbakhsh et al., 2015) and microbial consortia (Parkin and Murray, 2018). The former provides a particularly striking contrast, since the waves that drive Dictyostelium signaling are pulsatile in nature, yet are also used to coordinate chemotactic response. Whereas continuous emission relays create continuous, steep concentration profiles, pulsatile relays in Dictyostelium create traveling wave packets of high concentration, each of which elicits a chemotactic response. We see no evidence of ‘jumps’ in chemotactic response during neutrophil swarming. It is not clear what drives one organism to adopt pulsatile signaling over relays with continuous emission, or vice versa.

Finally, it would also be interesting to leverage the design principles we have discussed for engineering synthetic relays, a field with a rich history (Parkin and Murray, 2018; Brenner et al., 2008; Brenner et al., 2007; Basu et al., 2005). To that end, our results provide a general framework for determining how system dimensionality, diffusion constants, activation functions, cell density, etc. affect cell signaling and wave initiation. Experimental work on this problem and others would provide tests of our many quantitative predictions.

Materials and methods

To find the information wave front for cells in n dimensions and diffusion in m dimensions with continuous emission and Heaviside activation, we make use of the Green’s function for the diffusion equation with sources in n dimensions and diffusion in m dimensions, Gn,m(r,t;R,T). These equations are enumerated in Appendix 7: Initiation dynamics; dTdRaρGn,m(r,t;R,T) describes the concentration created at a radius r and time t by a tiny ring of sources at radius R with density ρ that emit at rate a for duration dT at time T.

To find the information front, one is looking for a curve rc(t) such that Cth=c(rc(t),t). Thus, with an initial signaling colony of size ri, one must solve the problem:

(13) Cth=aρ0t𝑑T0max[ri,rc(T)]𝑑RGn,m(rc(t),t;R,T).

This constraint equation considers every radius at time T and, if it is less than rc(T), adds a concentration contribution of aρdTdRGn,m(rc(t),t;R,T) at rc(t); the sum of all these contributions must be equal to Cth. If one wishes to find the information front for a simple diffusive theory, one performs the same integral as above, but truncates the integration over R at ri.

This method is preferable to brute PDE solving (for example, on a grid) since the former requires fine-grained meshing over the out-of-plane dimension when considering systems of, for example, cells in 2D and diffusion in 3D. In contrast, our Green’s function method requires only numerical integration over the in-plane sources; the Green’s functions appropriately keep track of the out-of-plane dynamics for us.

To solve this problem, we first find the initiation time, then find rc(t) at discrete times, incrementing in steps of ΔtD/v2 (we use Δt=D/10v2 in the main text and Appendices, which gives convergence of the information wave front). Linear interpolation between these points defines a continuous curve rc(t).

An explicit implementation of this method is provided at github./pdieterle/diffWavePropAndInit (Dieterle, 2020; copy archived at swh:1:rev:f8d9feffd57d05f47c8c14c6d9850643b2858d0a).

Appendix 1

Model set-up

Before doing any math, let’s set up the scenarios we intend to study. We consider a continuum of cells described by a cell density, ρ. These cells emit one type of signaling molecule at a rate a when the local concentration of the signaling molecule is above a certain threshold, Cth. The molecules diffuse in the extracellular medium with diffusion constant D. The concentration of the signaling molecule is described by the variable c, which is a function of both space and time: c=c(𝐫,t). In general, then, we have

(14) ct=D2c+aρΘ[c-Cth]

with Θ[.] the Heaviside step function. We study this model and variants going forward.

A word about notation before we proceed: in Equation (14), we have written the source term as aρΘ[c-Cth], with ρ the cell density. In some cases, we will write this source term slightly differently; for instance, with point-like cells homogeneously distributed in a two-dimensional plane, we will write the source term as aρδ(z)Θ[c-Cth]. Here, ρ is a two-dimensional cell density, δ(.) is the Dirac delta function, and z is the out-of-plane dimension. Therefore, the cell density ρ which appears in all of the subsequent discussion will be a three-dimensional density if the cells are distributed in three dimensions, a two-dimensional density if the cells are distributed in two dimensions, and an one-dimensional density if the cells are distributed in one dimension. As we will show below, the choice of a delta function to describe the distribution in out-of-plane dimensions is justified when the cell size is small compared to the inherent length scale of the diffusive wave, D/v.

Appendix 2

Asymptotic wave ansatz

We start out by seeking to understand what dynamical properties a system described by Equation (14) has at large times. To study these dynamics, we need an inspired guess for what the dynamics will look like. We imagine that when a small volume of cells starts signaling its neighbors – and those neighbors start signaling their neighbors – that a reasonable guess for the dynamics of such a signaling relay is an outward propagating wave with speed v. We define r as the distance from the center of the outward propagating wave. At long times for a uniform cell density, the shape information wave front will obey radial symmetry and our ansatz becomes c(𝐫,t)=c(𝐫-vtr^) with r^ the unit vector pointing from the origin to the wave front.

With this guess, we can define a new coordinate r~=r-vt (we call this x~=x-vt for cells in one dimension) which defines the distance to the wave front. Note that r~<0 means we are inside the wave front while r~>0 means we are beyond it. With these definitions, c/r~=c/r and c/t=-vc/r~. For cells in 1D, we consider the y and z to be dimensions perpendicular to the line of cells with the density described by ρδ(y)δ(z) with ρ measured in cells per unit length; for cells in 2D, we consider z to be the out-of-plane dimension and the density to be described by ρδ(z) with ρ measured in cells per unit area; for cells in 3D, ρ is measured in cells per unit volume. Assuming azimuthal symmetry in 2D and radial symmetry in 3D, we arrive at

(15a) cells in 1D: 0=D(2cx~2+2cy2+2cz2)+vcx~+aρ δ(y)δ(z)Θ[cCth]
(15b) cells in 2D: 0=D(2cr~2+1rcr~+2cz2)+vcr~+aρ δ(z)Θ[cCth]
(15c) cells in 3D: 0=D(2cr~2+2rcr~)+vcr~+aρ Θ[cCth]

These equations can be simplified once more by noting that we are considering asymptotic – that is, large r – dynamics. Thus, as long as vD/r, we can say that vc/r~ dominates terms like D(c/r~)/r and we can ignore the latter (Tanaka et al., 2017). By construction, our ansatz says that c(r~=0)Cth, meaning that Θ[c-Cth]=Θ[-r~]. This gives simplified equations according to:

(16a) cells in 1D: 0=D(2cx~2+2cy2+2cz2)+vcx~+aρ δ(y)δ(z)Θ[x~]
(16b) cells in 2D: 0=D(2cr~2+2cz2)+vcr~+aρ δ(z)Θ[r~]
(16c) cells in 3D: 0=D2cr~2+vcr~+aρ Θ[r~].

These equations provide both a natural length scale, D/v, and a natural timescale, D/v2. One can see that these are the relevant time and length scales for our problem by non-dimensionalizing, for example, Equation (16c) to get

(17) cells in 3D: 0=2(cv2/aρD)(vr~/D)2+(cv2/aρD)(vr~/D)+Θ[vr~/D].

Thus, every length scale in the problem is normalized by D/v; because of our traveling wave ansatz, every length scale can be converted to a time scale by dividing by v, giving D/v2 as the natural timescale. The natural length scale is useful, for example, for understanding what it means for cells to be in ‘one dimension’ or for diffusion to be in ‘two dimensions’. If the cells are organized in a line (or on a plane) such that their average deviation from the line (or distance from the plane) is dD/v, then they are effectively in one (or two) dimensions and Equation (16a) (or Equation (16b)) holds. If cells are constricted to a narrow channel of width hD/v (or an extracellular medium of thickness hD/v), then diffusion is effectively one (or two) dimensional. For instance, for cells confined in a very narrow one-dimensional channel of width hD/v, we can simplify Equation (16a) because 2c/y2=2c/z2=0 and δ(y)δ(z)1/h2. The resulting equation is the exact same as for cells in 3D but with a source term proportional to aρ/h2:

(18) cells in 1D, diffusion in 1D: 0=D2cx~2+vcx~+aρh2Θ[x~].

For cells in 2D with an extracellular medium of thickness hD/v, diffusion is effectively two-dimension and, by the same logic that produced Equation (18), the asymptotic governing equation is:

(19) cells in 2D, diffusion in 2D: 0=D2cr~2+vcr~+aρhΘ[r~].

As Equations (17), (18), and (19) are all the same, the dynamics of cells in 1D with diffusion in 1D are the same as those of cells in 2D with diffusion in 2D or those in 3D with diffusion in 3D.

Moreover, the non-dimensional Equation (17) shows us that the concentration scale of interest is aρD/v2, meaning that there must be a relationship along the lines of CthaρD/v2v(aρD/Cth)1/2 – exactly the wave speed relationship we found in the main text.

One can, however, arrive at a different governing equation by considering cells in 2D (for example, cells sitting on a plane) with a thick extracellular medium of thickness hD/v. In this case, diffusion effectively takes place in three dimensions and

(20) cells in 2D, diffusion in 3D: 0=D(2cr~2+2cz2)+vcr~+aρδ(z)Θ[r~].

which, by the same logic above, is functionally equivalent to the governing equation for cells in 1D with diffusion in 2D. We can therefore see that it is not the dimensionality of the cell distribution or the diffusive environment that determines the asymptotic dynamics, but rather the difference in dimension between the two.

Going forward, we will think of cells in two dimensions, as we have done in the main text. This will allow us to interpolate between an effectively two-dimensional diffusive environment (the thin extracellular medium limit) and an effectively three-dimensional environment (the thick extracellular medium limit).

Cells in 2D, diffusion in 2D: the thin extracellular medium limit

For an extracellular medium of thickness hD/v, diffusion effectively takes place in two dimensions as argued in the previous section. The signaling molecule concentration has no z-dependence and concentrations get normalized by h. Here,

(21) 0=D2cr~2+vcr~+aρhΘ[-r~]

is our asymptotic governing equation.

For both r~<0 and r~>0, Equation (21) reduces to two straightforward-to-solve linear ODEs. With bi as constants that we will determine momentarily,

(22a) r~<0: 0=D2cr~2+vcr~+aρ/h  c(r~<0)=b2evr~/D+b3aρr~/hv
(22b) r~>0: 0=D2cr~2+vcr~  c(r~>0)=b0evr~/D+b1.

We can solve for the bi by applying boundary conditions. First, we demand c0 as r~ and that the concentration profile only blow up linearly as r~-. Physically, these demands are justified as follows: the concentration as r~ has to go to zero because there are no cells emitting in that region and it is far from the wave front; the concentration as r~- can grow at most linearly because the cells a distance -r~ from the wave front have only been emitting for a time -r~/v. Combining these asymptotic boundary conditions with the demand that c(r~) be continuous and have a continuous first derivative at r~=0 allows us to stitch together the solutions in Equations (22a) and Equations (22b) to yield:

(23a) c(r~<0)=aρD/hv2-aρr~/hv
(23b) c(r~>0)=aρDe-vr~/D/hv2.

We show this concentration profile in the left panel of Appendix 3—figure 1A. From the above, we infer that

(24) Cth=c(0)=aρD/hv2  v=aρD/hCth.

Thus, we have an explicit formula relating wave speed, emission rate, cell density, diffusion constant, extracellular medium thickness, and threshold concentration. We can also see that the concentration profile beyond the wave front is exponential, not Gaussian as for simple diffusion. The concentration inside the wave front grows linearly as the distance from the wave front.

Cells in 2D, diffusion in 3D: the thick extracellular medium limit

Next, we consider Equation (16b) in the limit that the extracellular medium hD/v. In this limit, we effectively have cells in 2D with diffusion in 3D. With cells sitting on a substrate, signaling molecules can only diffusive in the upper half of the plane, and we have a semi-infinite environment which accounts for an extra factor of 2 in the emission term, yielding:

(25) 0=D(2cr~2+2cz2)+vcr~+2aρδ(z)Θ[-r~].

Instead of working directly with δ(z), we consider the cells to be of a thickness H such that 2δ(z)1H2πexp(-z2/2H2) and

(26) 0=D(2cr~2+2cz2)+vcr~+aρH2πe-z2/2H2Θ[-r~]

Next, we take a partial Fourier transform of Equation (26) with k and C(r~,k) the Fourier partners of z and c(r~,z), respectively. Here, we choose C(r~,k)12π-eikzc(r~,z). This gives:

(27) 0=D(2Cr~2-k2C)+vCr~+2πe-H2k2/2aρΘ[-r~]

which is another pair of piecewise, straightforward-to-solve, linear ODEs. Solving with the same boundary conditions that yielded Equations (23a) and Equations (23b), we arrive at

(28a) C(r~<0,k)=aρe-H2k2/2Dk22π(2-(1+4D2k2/v2+1)exp(vr~2D(4D2k2/v2+1-1))4D2k2/v2+1)
(28b) C(r~>0,k)=aρe-H2k2/2Dk22π4D2k2/v2+1-14D2k2/v2+1exp(-vr~2D(4D2k2/v2+1+1)).

To find the concentrations at the cells, we can take the inverse partial Fourier transform of these expressions at z=0. But first, we note that the right sides of Equations (28a) and Equations (28b) have no support when kv/D. Thus, if Hv/D1, the term e-H2k2/2 is irrelevant for calculating the real-space concentrations and can be replaced with 1. This is equivalent to having chosen δ(z) to describe the out-of-plane cell density.

Proceeding with e-H2k2/21, one can take the inverse partial Fourier transform and arrive at

(29) Cth=2aρ/πv.

Similarly, in the limit |r~|D/v,

(30a) c(r~-D/v,z)2aρv-r~vπD(evz2/4Dr~--πvz24Dr~erfc-vz24Dr~)
(30b) c(r~D/v,z)aρDπr~v3e-vr~/De-vz2/4Dr~.

The former has the same functional dependence on z as the concentration a distance z away from a continuously emitting point source with diffusion in 1D after a time r~/v (see Appendix 6: Assessing the validity of a continuum analysis). Using the above, we can find the concentration in the plane of the cells (z=0):

(31a) c(r~D/v,z=0)2aρv-r~vπD
(31b) c(r~D/v,z=0)aρDπr~v3e-vr~/D.

We show this concentration profile in the left panel of Appendix 3—figure 1B.

Of course, one can arrive at the scaling form of Equation (29) through dimensional analysis considerations, as discussed in the main text. Equivalently, one can non-dimensionalize Equation (25) along the lines of Equation (17), as we will now. By normalizing all length scales by D/v and all concentration scales by aρ/v, one can show that Equation (25) is equivalent to

(32) 0=(2(cv/aρ)(r~v/D)2+2(cv/aρ)(zv/D)2)+(cv/aρ)(r~v/D)+2δ(zv/D)Θ[-vr~/D].

from which we can tell that the concentration scales of the problem, including Cth, are proportional to aρ/v, and thus that vaρ/Cth, exactly as required by Equation (29).

Cells in 1D, diffusion in 3D: an artificial case

Finally, we consider a line of cells in one dimension with diffusion taking place in three dimensions. This corresponds to a somewhat artificial test case of cells in a line with mean distance from the line dD/v and diffusion in an environment of size hD/v in the dimensions perpendicular to this line of cells. Nonetheless, it is interesting because we have to include the finite size of the cells in order to get a traveling wave solution. Here,

(33) 0=D(2cx~2+2cy2+2cz2)+vcx~+aρπH2e-(y2+z2)/2H2Θ[-r~]

with v the wave speed; x~ the distance to the wave front; y and z the extra diffusive dimensions; and H the size of the cells.

Taking partial Fourier transforms across both y and z gives with the Fourier transform conventions and notation used above gives

(34) 0=D(2Cx~2-ky2C-kz2C)+vCx~+aρπe-(ky2+kz2)H2/2Θ[-x~]

which reduces to Equation (27) with k2ky2+kz2. One can then find the concentration profiles and self-consistency relationship for Cth by inverse Fourier transforming the analogs of Equations (31a) and Equations (31b). The value of Cth=c(x~=0,y=z=0) diverges as H0.

Appendix 3

Pulsed emission and decay

In this section, we consider pulsed emission and decay of the signaling molecule. These scenarios are relevant for signaling pathways in, for example, Dictyostelium (Pálsson and Cox, 1996; Noorbakhsh et al., 2015; Kessler and Levine, 1993) and E. coli (Parkin and Murray, 2018), in which intracellular dynamics produce a pulse-like release of signaling molecules into the extracellular medium.Kessler and Levine, 1993 have previously used this machinery to construct a signaling model for Dictyostelium, including pulsed emission and signaling molecule decay. Here, we consider the effects of each independently.

We explicitly discuss only the asymptotics of cells in 2D with diffusion in 2D (equivalent to cells in 1D with diffusion in 1D or cells in 3D with diffusion in 3D, as shown previously), although we quote the results for cells in 2D with diffusion in 3D (equivalent to cells in 1D with diffusion in 2D) which are obtained using the Fourier transform machinery in Appendix 2: Asymptotic wave ansatz. Here again, the asymptotic dynamics depend on the difference in dimensionality between the cellular and the diffusive environment.

Pulsed emission with cells in 2D and diffusion in 2D

Here, we consider a square pulse of length τ emitted once a cell exceeds the threshold concentration Cth. In the moving frame, this pulse has length vτ – for notational simplicity here, we dispense with dimensional subscripts on the wave speed – giving rise to the pulsed emission analog of Equation (21):

(35) 0=D2cr~2+vcr~+aρhΘ[-r~]Θ[r~+vτ]

which is nothing more than three piecewise linear equations, which we stitch together as before. The source term is zero when r~<-vτ (Region I) or r~>0 (Region III) and aρ0 for -vτ<r~<0 (Region II). We thus recover the following:

(36a) Region I: c(r~<-vτ)=b1+b2e-vr~/D
(36b) Region II: c(-vτ<r~<0)=b3+b4e-vr~/D-aρr~/hv
(36c) Region III: c(r~>0)=b5+b6e-vr~/D

Applying the same boundary conditions as with continuous emission, we arrive at

(37a) Region I: c(r~<-vτ)=aρτ/h
(37b) Region II: c(vτ<r~<0)=aρD/hv2aρDhv2ev2τ/Devr~/Daρr~/hv
(37c) Region III: c(r~>0)=aρDhv2(1-e-v2τ/D)e-vr~/D

which tells us a few things of interest. First – as seen in the right panel of Appendix 3—figure 1 A (here, τ=2D/v2=2Cth/aρ) – the concentration profile for r~<-vτ is flat. (As with continuous emission, pulsed emission gives the familiar exponential profile beyond the wave front.) Second, the wave speed, pulse width, cell density, emission rate, extracellular medium thickness, and threshold concentration are related through the equation

(38) Cth=aρD(1-e-v2τ/D)/hv2.
Appendix 3—figure 1
Relay concentration profiles.

(A) Signaling molecule concentration profiles for cells in 2D and diffusion in 2D. Here, we assume continuous emission and no signaling molecule decay (left panel), continuous emission and signaling molecule decay (middle panel), or pulsed emission with no signaling molecule decay (right panel). The decay constant for the middle panel is γ=v2/4D while the pulse width is τ=2D/v in the right panel. For the left, center, and right panels, the threshold concentration is calculated according to Equations (24), (41), and (38), respectively. As discussed in the main text, the concentration profiles flatten (with respect to the profile generated by continuous emission without decay) inside the wave front once decay or pulsed emission is accounted for. (B) Signaling molecule concentration profiles for the same cases as in A, but with diffusion in 3D. Compared to the case of continuous emission without decay, the concentration profiles flatten (when accounting for decay) or have a local maximum (in the case of pulsed emission).

For τD/v2, we recover the usual relationship of Cth=aρD/hv2. In region 2, the profile will grow linearly as before until -vr~/D becomes comparable to v2τ/D, at which point e-v2τ/D-vr~/D becomes of order unity and the profile levels off.

To understand how the wave speed with pulsed emission, v, compares to the wave speed with continuous emission, (aρD/hCth)1/2, we have plotted v/(aρD/hCth)1/2 as a function of 1/τ in Appendix 3—figure 2A. We have normalized τ by a characteristic time τc=hCth/aρ, which is equal to D/v2 for continuous emission. When τ<τc, the wave speed goes to zero. There is no wave-like solution for shorter pulses.

Appendix 3—figure 2
Wave speed with pulsed emission or decay.

(A) Wave speed v for square pulse emission by cells in 2D with diffusion in 2D as a function of pulse width τ. We normalize v by the τ wave speed of (aρD/Cth)1/2. At τ=τc=hCth/aρ, the wave speed goes to zero, and for shorter pulses there is no wave-like solution. (B) Wave speed for continuous emission by cells in 2D with diffusion in 2D, accounting for signaling molecule decay at rate γ. For γ=tmin, 2D=1/2τc, with τc as in panel A, the wave speed goes to zero and there are no wave-like solutions for larger decay rates. (C) Wave speed v as a function of pulse width τ for square pulse emission of a signaling molecule by cells in 2D with a 3D diffusive environment. At τ1.53τc=1.53D(2Cth/πaρ)2 (vertical dashed line), there is a minimum wave speed of v0.6×(2aρ/πCth) (horizontal dashed line), unlike with 2D diffusion. Importantly, τ1.53τc is longer than the minimum initiation time of 4τc/π. Thus, for values of τ between these two, cells can reach the threshold concentration but cannot propagate a traveling information wave. (D) Wave speed v as a function of decay rate γ for cells in 2D with diffusion in 3D. At γ=(π/4)2/τc, with τc as in C, the wave speed goes to zero.

We note that a timed pulsed emission considered here is formally equivalent to cells signaling until the local concentration exceeds c(r~=-vτ)=aρτ/h. This is relevant in, for example, quorum sensing models in which the local presence of a signaling molecule can both upregulate (at relatively low concentrations) and downregulate (at relatively high concentrations) release of the same signaling molecule (Parkin and Murray, 2018).

Continuous emission plus decay with cells in 2D and diffusion in 2D

At last, we characterize the effect of signaling molecule decay at rate γ by adding a term of -γc to Equation (16b). In the thin extracellular medium limit,

(39) 0=D2cr~2+vcr~+aρhΘ[-r~]-γc.

This is another piecewise set of linear differential equations, which we can solve as without decay to yield:

(40a) c(r~<0)=aρhγ-aρ2hγ[1+(1+4Dγ/v2)-1/2]er~v2D(4Dγ/v2+1-1)
(40b) c(r~>0)=aρ2hγ[1-(1+4Dγ/v2)-1/2]e-r~v2D(4Dγ/v2+1+1)

as the concentration profiles and

(41) Cth=aρ2hγ[1-(1+4Dγ/v2)-1/2]v=4Dγ[(1-2hγCthaρ)-2-1]-1

as our wave speed relationship. The concentration profile is flatter than its decay-free counterpart (Appendix 3—figure 1A). For γv2/D, Equation (41) gives the decay-free relationship. And – as with pulsed emission – the wave speed goes to zero, this time when γ1/2τc where τc=hCth/aρ (Appendix 3—figure 2B).

Pulsed emission with cells in 2D and diffusion in 2D


(42) 0=D(2Cr~2-k2C)+vCr~+2/πaρΘ[-r~]Θ[vτ+r~]

which we can solve to yield:

(43a) Region I: C(r~<-vτ,k)=B1(k)er~v2D(-1+1+4D2k2/v2)
(43b) Region II: C(vτ<r~<0,k)=2πaρ0Dk2+B2(k)er~v2D(1+1+4D2k2/v2)+B3(k)er~v2D(1+1+4D2k2/v2)
(43c) Region III: C(r~>0,k)=B4(k)e-r~v2D(1+1+4D2k2/v2)

with Bi(k) chosen such that C and its first derivative are continuous:

(44a) B1(k)=aρD2π(v+4D2k2+v2)(ev2τ2D(4D2k2/v2+1-1)-1)k24D2k2+v2
(44b) B2(k)=aρD2π(v-4D2k2+v2)e-v2τ2D(4D2k2/v2+1+1)k24D2k2+v2
(44c) B3(k)=-aρD2πv+4D2k2+v2k24D2k2+v2
(44d) B4(k)=aρD2π(v-4D2k2+v2)(1-e-v2τ2D(4D2k2/v2+1-1))k24D2k2+v2

Inverse Fourier transforming at z=0 gives the real-space concentration and the following self-consistency relationship for the wave speed:

(45) Cth=12π-𝑑kB4(k),

which simplifies to Cth=2aρ/πv for τD/v2.

We can compare τ to the characteristic time τc=D(πCth/2aρ)2, which is equal to D/v2 in the limit of continuous emission. Numerical solution of Equation (45) reveals that there is no self-consistent solution to Equation (45) until about τ1.53τc, at which point v0.6×2aρ/πCth;τ1.53τc is larger than the minimum initiation time of tmin,3D=4τc/π (Appendix 3—figure 2 C, Section Appendix 7: Initiation dynamics). Thus, an initial pulse of length 4τc/π<τ<1.53τc from cells within an initial signaling radius ri can cause neighboring cells to exceed Cth, but cannot trigger a wave-like solution asymptotically.

Continuous emission plus decay with cells in 2D and diffusion in 3D

We can add signaling molecule decay to the embedded system dynamics by adding a term of -γC(r~,k) to Equation (27) with γ the signaling molecule decay rate. Going through the same exercise yields:

(46a) C(r~<0,k)=2πaργ+k2D-aρ2π4D2k2/v2+4Dγ/v2+1+1(Dk2+γ)4D2k2/v2+4Dγ/v2+1er~v2D(-1+4D2k2/v2+4Dγ/v2+1)
(46b) C(r~>0,k)=aρ2π4D2k2/v2+4Dγ/v2+1-1(Dk2+γ)4D2k2/v2+4Dγ/v2+1e-r~v2D(1+4D2k2/v2+4Dγ/v2+1)


(47) Cth=aρπDγarcsin[(1+v2/4Dγ)-1/2]

as the parameter relationship obtained after an inverse Fourier transform at z=0 and x=0. This gives a profile that propagates as a pulse (Appendix 3—figure 1B).

In the limit γv2/D, we recover the familiar expression Cth=2aρ/πv. Again, as when we accounted for signaling molecule decay with cells in 2D and diffusion in 2D, the wave speed approaches zero, but with γ=(π/4)2/τc where τc=D(πCth/2aρ)2 (Appendix 3—figure 2 D).

Appendix 4

Finite extracellular medium

We consider now what happens when one does not lie in the extreme cases of an extracellular medium of thickness hD/v or hD/v. To examine the case of arbitrary thickness h, we turn to the method of images (Appendix 4—figure 1A).

Appendix 4—figure 1
Wave speed in a finite extracellular medium.

(A) Method of images configuration for solving for concentrations in a finite-sized extracellular medium of thickness h. A point-like source (a cell) at z=0 is placed in an extracellular medium of thickness h. In order to calculate the concentration with the appropriate boundary conditions of c/z at z=0,h, one need only add the contributions from ‘image cells’ at z=±2jh for j=1,2,. (B) Universal curve (black line) showing non-dimensionalized wave speed (vh/2D) versus non-dimensionalized threshold concentration (πCthD/aρh). In the limit of vh/2D1, we recover the familiar 3D scaling law of Equation (29) (blue line). In the limit of vh/2D1, we get the 2D scaling law of Equation (24) (red line).

Our boundary conditions for the extracellular medium require the concentration to obey c/z=0 at z=0,h – signaling molecules cannot diffuse through the boundaries. By invoking the uniqueness theorem, we know that if we can find an arrangement of 'image cells' – each emitting diffusible signaling molecule at rate a – that satisfies these boundary conditions, then this arrangement of cells gives the unique solution for the concentration profile inside the extracellular medium. In our case, to satisfy the boundary condition above, we have image cells at z=±2jh for j=1,2,.

This means that we can find the concentration profiles simply by adding up the contributions from many discrete sources. Given this knowledge, we seek a relationship like Equation (24) or Equation (29) but for arbitrary h. To do so, we use Green’s function integration and the fact that we can analyze the asymptotic dynamics of cells in 1D to deduce the asymptotic dynamics of cells in 2D, as previously shown.

The concentration – as measured at (r,z=0,t) – of a burst-like emission by a single point-like source at (R,2jh,T) is given by the Green’s function:

(48) G(r,z=0,t;R,2jh,T)=e-(r-R)2+(2jh)24D(t-T)2πD(t-T),

so the Green’s function of a single point-like source in a finite-thickness extracellular medium is (Appendix 4—figure 1A):

(49) Gh(r,z=0,t;R,T)=j=-G(r,z=0,t;R,2jh,T)=e-(r-R)24D(t-T)2πD(t-T)(1+2j=1e-(jh)2D(t-T)).

We assume a traveling wave solution at speed v meaning that the concentration Cth=c(r=vt,z=0,t) is, for a density of cells ρ emitting with rate a, given by:

(50) Cth=c(vt,0,t)=aρlimt0vtdRR/vtdTGh(vt,0,t;R,T)=aρ2πD-0𝑑R~0-R~/vdt~t~e-R~24Dt~(1+2j=1e-(jh)2Dt~)

with the substitutions t~=t-T and R~=R-vt. This yields:

(51) Cth=2aρπv(1+2j=10-𝑑xEi[1x(jhv2D)2+x])=2aρπvj=-0-𝑑xEi[1x(jhv2D)2+x]

for x=vR~/4D and Ei[.] the exponential integral function. For hD/v, Equation (51) reduces to Equation (29) because the term j=10. Meanwhile, for hD/v, the sum over j can be turned into an integral, giving the familiar thin extracellular medium relationship, Equation (24).

We emphasize that Equation (51) provides a universal relationship between threshold concentration, wave speed, cell density, and signaling molecule emission rate for any extracellular medium thickness. By dividing both sides of Equation (51) by aρh/πD, we arrive at a relationship between a non-dimensionalized threshold concentration, πCthD/aρh, and a non-dimensionalized wave speed, vh/2D:

(52) πCthDaρh=2Dvhj=-0-𝑑xEi[1x(jhv2D)2+x].

We plot this relationship in Appendix 4—figure 1B and see that Equation (52) is an interpolation between the thin (hD/v) and thick (hD/v) extracellular medium limits.

Appendix 5

Asymptotic wave dynamics with Hill function activation

Numerical solutions show traveling waves

As shown in the main text, making the change from a Heaviside function source term to an order-n Hill function (Θ[c-Cth]cn/(cn+Cthn)) in Equation (14) preserves the scaling relationships Equations (24) and (29) with a constant factor as long as the new source terms give traveling wave solutions. We have found numerically that n1 Hill functions indeed give traveling wave solutions, with n=1,2,3 shown for thin and thick extracellular media in Appendix 5—figure 1.

Appendix 5—figure 1
Wave speeds and profiles with Hill function activation.

(A) Numerical simulation of cells in 1D with diffusion in 1D and Hill function activation. The details of the simulation are described in Appendix 5: Asymptotic wave dynamics with Hill function activation. For n2, we see good agreement between the Heaviside theory (black lines) and the Hill function numerics (red dots). Snapshots are shown at t/τc=5,10,15,20,25 where τc=h2Cth/aρ equals D/vΘ2, the characteristic time scale for Heaviside activation. Note that the x-axis is scaled by the characteristic length lc=h(DCth/aρ)1/2, which is the length scale for Heaviside activation with a delta function source, D/vΘ. In the insets, we display the wave speed for the order-n Hill function, vn, compared to the Heaviside activation wave speed, vΘ=(aρD/h2Cth)1/2. Our fit to the n=1 data gives vn=11.93vΘ, but we have shown in Appendix 5: Asymptotic wave dynamics with Hill function activation that vn=1=2vΘ, meaning that the wave speeds in the insets are slight underestimates. (B) Numerical simulation of cells in 1D and diffusion in 2D with Hill function activation. For n2, the wave speed and concentration profiles (blue dots) agree well with the Heaviside theory (black lines). The theory plotted here assumes a delta-function-like source with respect to the extra diffusive dimension, as in Equation (25). The numerics, however, use a very narrow (H=lc/10) Gaussian source. Here, the characteristic length lc=πhDCth/2aρ equals D/vΘ, the length scale for Heaviside activation. The x-axis is scaled by the same quantity. Snapshots are shown at t/τc=5,10,15,20,25 with the characteristic time τc=D(πhCth/2aρ)2 equal to D/vΘ2, the time scale for Heaviside activation. In the insets, we display the wave speed for the order-n Hill function, vn, compared to the Heaviside activation wave speed, vΘ=2aρ/πCth.

To find these solutions, we numerically solved Equation (14) with a Hill function source term for cells in 1D with diffusion in one or two dimensions. We used D=10-10 m2/s and vΘ=2 µm/s with the threshold concentration determined by Equations (29) and (24) with vvΘ. In this way, we could compare vn – the wave speed given by the order-n Hill function – to the Heaviside wave speed vΘ. For our numerics, we imposed a maximum size step of D/10vΘ for the spatial dimension where the cells live, a maximum step size of D/30vΘ for the added spatial dimension (when modeling diffusion in 2D), and a maximum time step of D/10vΘ2. These step sizes give convergence of the information wave fronts, which we define as the curves rc(t) such that c(rc(t),t)=Cth. We simulated times tmax25D/vΘ2 using Mathematica’s ‘NDSolveValue’ function and, when modeling diffusion in 2D, replaced the delta function in Equation (14) with a Gaussian of width D/10vΘ. As noted earlier, since the width of this Gaussian is sufficiently smaller than D/v, the numerical results will be close to the delta function limit; indeed, this substitution gives a wave speed (according to an inverse Fourier transform of Equation (28b) at r~,z=0) of v/vΘ0.95, in close correspondence with the Delta function wave speed. The initiating colony is of size ri=2D/vΘ, and we assume all cells in the initiating colony signal at the maximal rate a.

To find the wave speeds noted in Equation (1), we found the location of the wave front, rc(t), and fit a line to the region of 24D/vΘ2t25D/vΘ2. As shown in Appendix 5—figure 1, for n2, the wave speeds are very close to vΘ with significant deviation only for n=1. Even the concentration profiles are in good agreement with the Heaviside solution for n2.

Connection to Fisher waves

When the dimensionality of the cell distribution matches the diffusive dimensionality and cell activation is described by the n=1 Hill function, one can find the wave speed by using a modified version of the analysis pioneered by Fisher, 1937 and Kolmogorov et al., 1937. To see this, we first consider a modified version of Equation (16c) with order-n Hill function activation:

(53) 0=D2cr~2+vcr~+aρcncn+Cthn.

The new activation term is mathematically obnoxious as it no longer has a simple spatial interpretation; with a Heaviside function, for example, one can turn a term like Θ[c-Cth] into a simple function of r~:Θ[cCth]=Θ[r~]. This has an important consequence: instead of solving two differential equations with constant source terms and matching boundary conditions (as we did for the Heaviside emission), we must now solve a single differential equation with a difficult non-linear source term. We note that Equation (53) is, in the limit n, equivalent to a relay with Heaviside activation.

However, there is a distinct advantage to the new source term: it is one-to-one in c. Thus, we may make the substitution f=cncn+Cthnc=Cth(f1-f)1/n, then plug this into Equation (53) to yield (after some rearrangement):

(54) 0=vfr~+D[2fr~2+1+(2f-1)nnf(1-f)(fr~)2]+naρCthf2-1/n(1-f)1+1/n.

Equation (54) looks a lot like the traditional Fisher equation,

(55) 0=vfr~+D2fr~2+naρCthf(1-f)

in that it has a source term that goes to zero at f{0,1}, a term D2f2r~, and a term vfr~. We therefore take Fisher’s approach (Fisher, 1937) and think of gradients of f as functions of f rather than x. As such, we define F(f)=fx, which allows us to make the substitution 2fx2=Fx=fxFf=FFf. This is valid under the assumption that the concentration profiles are monotone decreasing, which is seen to be the case in Appendix 5—figure 1. Under all of the above, we get:

(56) 0=vF+D[FFf+1+(2f-1)nnf(1-f)F2]+naρCthf2-1/n(1-f)1+1/n

which gives us a non-linear ODE for F. Of particular note are the boundary conditions for F. Namely, F(0)=F(1)=0, which is to say that cells well inside of the wave front are all emitting at their maximal rate since cCth and that cells well beyond the wave front are not emitting at all. Thus, there is no spatial dependence on the cellular activation, f=cn/(cn+Cthn), in these regions.

Next, we turn to the traditional method of examining the f0 limit. As F0, Equation (56) becomes, to lowest order in f and assuming F-λfβ,

(57) 0=-λv+D[λ2βfβ-1+1-nnλ2fβ-1]+naρCthf2-1/n-β.

where λ>0. One can only obtain a self-consistency relationship between v and Cth/aρ if β=2-1/n. Otherwise, f2-1/n-β diverges or goes to zero. With this choice, Equation (57) becomes:

(58a) n=1: 0=-λv+Dλ2/n+aρCthλ=12D(v±v2-4aρDCth)
(58b) n>1: 0=-λv+naρ0Cth+𝒪(f1-1/n)λ=naρvCth

where in Equation (58b) we have taken the f0 limit. To get a wave speed, v, out of Equation (58a), we demand that the quantity under the square root be non-negative, which ensures that λ>0 is a real number as assumed. This means v2aρD/Cth=2vΘ – a bound that is very similar to conventional Fisher waves. In the case of Fisher waves, the minimum wave speed is selected for (Fisher, 1937; Kolmogorov et al., 1937); the same is true here, as the minimum wave speed v=2aρD/Cth is what one finds after numerically solving the 1D dynamics (Appendix 5—figure 1A). In contrast, for n>1, this method yields no such wave speed bound.

Appendix 6

Assessing the validity of a continuum analysis

Next, we consider the validity of a continuum analysis like Equation (16b) for studying asymptotic wave dynamics. To do so, we compare our continuum wave speed relationships, Equations (24) and (29), to a simple model of discrete cells on a lattice in 1D (Appendix 6—figure 1A). We refer to this as the discrete lattice model, which has been studied previously and is the subject of on-going work (Keener, 2000; Dieterle P and Amir A, 2020. Manuscript in preparation). We briefly discuss the results of this lattice model below, and show that it agrees with the continuous model when the separation between cells, d, is much less than the characteristic length D/v. This heuristic also holds for cells in two and three dimensions, and for cells scattered randomly according to a Poisson process.

Appendix 6—figure 1
Wave dynamics with discrete sources.

(A) Cartoon showing the discrete lattice model of cell signaling relays for cells in 1D. Here, a concentration wave propagates rightward at speed v through a group of cells spaced by a constant distance d. The concentration at x~=t=0 defines the threshold concentration in this discrete system and is a function of v,D, and d: Cth=Cth(v,D,d). (B) Wave speed v compared to the continuous theory value of (aD/dh2Cth)1/2 as a function of vd/4D for cells in 1D with diffusion in 1D. We can see that the wave speed approaches the continuous theory value for vd/4D1. (C) Same as B but for cells in 1D and diffusion in 2D.

To start, we first calculate the three-dimensional concentration (SI units of 1/m3) generated by a continuously emitting point source emitting at a rate a at a distance x after a time t. For diffusion in m dimensions, we will refer to this concentration as c,m(x,t). For diffusion in m=2 dimensions, we will consider a semi-infinite environment in order to recapitulate Equation (29), which holds for cells in one (two) dimensions with diffusion in a semi-infinite two-dimensional (three-dimensional) space. These relationships are:

(59a) c,1(x,t)=ah20t𝑑Te-x2/4DT(4πDT)1/2=ah2tπD(e-x2/4Dt-πx24Dterfcx24Dt)
(59b) c,2(x,t)=2ah0t𝑑Te-x2/4DT4πDT=-a2πhDEi(-x24Dt)

with erfc[.] the complementary error function and Ei[.] the exponential integral. We have assumed one-dimensional diffusion takes place in a narrow channel with cross-sectional area h2 and two-dimensional diffusion takes place in an extracellular medium of thickness h.

Next, we assume the cells in this lattice model perform a signaling relay with Heaviside activation: once the local concentration exceeds a threshold Cth, they participate in the signaling molecule emission at a rate a. If the resulting wave speed is v, then a cell at a distance x~=-jd from the wave front has been emitting for a time jd/v. That cell then creates a concentration cj,m=c,m(jd,jd/v) at (x~,t)=(0,0). The full concentration at the wave front is Cth by definition, and it is equal to the sum of the concentrations created by all cells behind the wave front. This gives us the following self-consistency relationships between the threshold concentration, Cth; wave speed, v; diffusion constant, D; and cell separation d:

(60a) diff. in 1D: Cth=j=1cj,1=avh2j=1jvdπD[e-jvd/4D-jπvd4Derfc(jvd4D)]
(60b) diff. in 2D: Cth=j=1cj,2=-a2πhDj=1Ei(-jvd4D).

Equations (60a) and (60b) provide relationships analogous to Equations (24) and (29). In fact, in the limit vd/4D1, the sums in these relationships are well-approximated by an integral over j from j=0 to . In this limit, with ρ=1/d, Equation (60a) becomes Cth=aρD/h2v2, the one-dimensional analog of Equation (24); similarly, Equation (60b) simplifies to Cth=2aρ/πhv, the one-dimensional analog of Equation (29). (One can turn Equation (60b) into an integral from j=0 to . The integrand diverges at j=0, but is still integrable because the divergence is logarithmic.) Therefore, we can see that the continuum limit is valid when the separation between cells is d4D/v. We have shown the approach to the continuous theory limit in Appendix 6—figure 1B/C.

Appendix 7

Initiation dynamics

In this section, we demonstrate the initiation time relationships discussed in the main text using Green’s function integration. To do so, we write down the Green’s functions Gn,m(r,t;R,T) describing diffusion of molecules in m dimensions released by cells in n dimensions at (R,T) and measured by cells at (r,t). For nm, we assume a semi-infinite environment. For cells in 1D and diffusion in 1D, we assume a narrow channel of width h in both dimensions perpendicular to the channel. For cells in 1D or 2D and diffusion in 2D, we assume an extracellular medium of thickness h. We calculate the Green’s functions for cells in two dimensions by integrating over a ring of diffusive sources at radius R; we calculate the Green’s functions for cells in three dimensions by integrating over a shell of diffusive sources at radius R. Below, I0[.] is the zeroth I-Bessel function and sinh [.] is the hyperbolic sine function. The Green’s functions are

(61a) G1,1(r,t;R,T)=e-(r-R)2/4D(t-T)/h24πD(t-T)
(61b) G1,2(r,t;R,T)=e-(r-R)2/4D(t-T)/2πhD(t-T)
(61c) G2,2(r,t;R,T)=RI0[rR/2D(t-T)]e-(r2+R2)/4D(t-T)/2hD(t-T)
(61d) G2,3(r,t;R,T)=RI0[rR/2D(t-T)]e-(r2+R2)/4D(t-T)/2πD3/2(t-T)3/2
(61e) G3,3(r,t;R,T)=Rsinh[rR/2D(t-T)]e-(r2+R2)/4D(t-T)/r[πD(t-T)]1/2.

One-by-one, we study the initiation time for these systems by studying the self-consistency relationship for the threshold concentration Cth and initiation time tinit for a given initial signaling colony of size ri:

(62) Cth=aρ0tinit𝑑T0ri𝑑RGn,m(ri,tinit;R,T)=aρ0tinit𝑑T0ri𝑑RGn,m(ri,0;R,-T)

where the logic here is that the signaling wave initiates when the threshold concentration at the edge of the initial signaling colony exceeds Cth. At tinit, cells outside the colony participate in the signaling and birth a diffusive wave with dynamics we have already studied extensively. This scenario assumes the cells do not move – that the cell density is fixed. In the case of neutrophils, this means that we are ignoring the possibility that a cell initially located off the target randomly encounters the target and starts signaling. Unlike the asymptotic dynamics, the difference in diffusive and cell dimension is not the salient parameter for understanding diffusive wave initiation. Rather, the initiation dynamics are determined solely by the dimension of the diffusive environment.

We now seek to derive the equations in the main text which show the relationship between the wave initiation time tinit and initial signaling colony size ri for various system dimensionalities when riD/v or riD/v. The full relationships of tinit versus ri for various system dimensionalities are plotted in Figure 3 of the main text.

Initiation with cells in 1D and diffusion in 1D

With cells and diffusion in 1D, we can perform the integral Equation (62) directly (for cells in 1D, we consider the bounds on the integral over R to be -ri to ri) and get a closed form relationship:

(63) Cth=aρri2h2D[(Dti/πri2)1/2e-ri2/Dti-1+(1+Dti/2ri2)erf(ri2/Dti)1/2].

In the limit where riDti, we get that ti=2D/v2 by using the asymptotic relationship Equation (24). This is the minimum initiation time, tmin,1,1.

(64) tmin,1,1=2D/v2

and it tells us that riDti is equivalent to riD/v – we can appeal to the natural length and time scales from our asymptotic analysis. We will soon see that this is also the minimum initiation time for cells in 2D with diffusion in 2D and cells in 3D with diffusion in 3D; this is the case because, as in the asymptotic analysis, we’ve essentially ignored the curvature of the target when calculating the riD/v initiation time.

In the opposite limit – riD/v – we can Taylor Expand Equation (63) and get

(65) riD/v:tinit(πD/4v2)(D/vri)2,

thus validating our equations in the main text.

Initiation with cells in 1D and diffusion in 2D

Next, we consider the self-consistency equation Equation (62) with n=1,m=2. As before, we first consider the limit of riDv and recover (through Equation (29)):

(66) tmin,1,2=4D/πv2

while for riDv,

(67) riD/v:log(Dtinit/ri2)2D/vri.

Initiation with cells in 2D and diffusion in 2D

Moving on, we consider the case in the main text of cells in 2D with diffusion in 2D. To perform the integration of Equation (62) in this case, it is easiest to rewrite the Bessel function in Equation (61c) in integral form, then integrate first over time. With riD/v, such an analysis gives a minimum initiation time of

(68) tmin,2,2=2D/v2.

In the opposite limit of riD/v,

(69) riD/v:log(4Dtinit/ri2)(2D/vri)2

as noted in the main text.

Initiation with cells in 2D and diffusion in 3D

Now, we will see that diffusive waves do not always initiate in 3D environments. We consider the integral in Equation (62), but take tinit which gives us a maximum concentration Cmax,2,3 at r=ri of:

(70) Cmax,2,3=2aρri/πD

Thus, by Equation (29), if ri<D/v, Cmax,2,3<Cth and the signaling wave cannot initiate. Examining Equation (62) for riD/v reveals

(71) riD/v:tinit(πD/16v2)(vri/D)4(vri/D-1)-2


(72) tmin,2,3=4D/πv2.

Initiation with cells in 3D and diffusion in 3D

Finally, we consider initiation with cells in 3D. Again, we consider the limit tinit in Equation (62) to get:

(73) Cmax,3,3=aρri2/3D.

Thus, waves do not initiate below a critical ri. However, here, the critical value is ri=3D/v. As with 1D cells/diffusion and 2D cells/diffusion, we recover

(74) tmin,3,3=2D/v2

in the limit riD/v.

Wave initiation with discrete cells

In the main text, we claimed that the qualitative wave initiation findings hold even in systems with discrete cells. Here, we show this by considering the extreme case of a single, point-like cell at the center of an initiating colony of size ri. We consider cells and diffusion in the same number of dimensions, m. In an m-dimensional diffusive environment, the concentration created by this source – which emits at a rate a – at a radius ri and time t will be given by:

(75) c(ri,t)=a0t𝑑Te-ri2/4DT(4πDT)-m/2.

For m=3, this integral is bounded from above by a/4πDri as t while for m=1,2, the concentration diverges as t. Thus we can see that, as in the continuum theory of the main text, the wave will always initiate for m=1,2 but has a critical radius for m=3.

However, this is not the only qualitative similarity between the discrete and continuum cases. For m=1,2, we obtain

(76a) cm=1(ri,t)=a4πD[2te-ri2/4Dt-riπ/Derfcri4Dt]
(76b) cm=2(ri,t)=a4πDΓ0(ri2/4Dt),

which we can manipulate to a more familiar form by considering the riD/v and tD/v2 limit, in which case we yield:

(77a) riD/v:cm=1(ri,t)atπD
(77b) riD/v:cm=2(ri,t)a4πDlog4Dtri2.

Making the substitution of c(ri,t)=Cth, setting t=tinit, multiplying both sides of the equations by rim, and writing Cth=aρD/v2aD/rim, we arrive at:

(78a) riD/v,m=1:tinit(D/v2)(D/vri)2
(78b) riD/v,m=2:log(Dtinit/ri2)(Dri/v)2,

exactly the results of Equations (65) and (69).

Thus, not only the rough phenomenology (of m=3 being distinct from m=1,2 in that exhibits a critical radius), but also the specific scalings of initiation time hold when accounting for the discreteness of cells.

Appendix 8

Wave initiation with Hill function activation

We now characterize the effect of Hill function-like activation on the signaling wave initiation, as we have done already for asymptotic dynamics. To do so, we consider a simple situation: cells within a volume of size ri signal with some rate a while neighboring cells outside of the initial signaling volume (i.e., with r>ri) signal with a concentration-dependent rate of acn/(cn+Cthn). These calculations give us an idea of how sensitive the initiation dynamics are to the details of cell activation. As we will show, the initiation dynamics with Hill function activation are a good approximation of the initiation dynamics for Heaviside activation when for relatively small n. One can imagine that such a situation may be relevant in, for example, the neutrophil swarming experiments presented in the main text (Reátegui et al., 2017). Here, cells in direct contact with a foreign protein begin signaling their neighbors, which respond to the presence of the signaling molecule by participating in the emission themselves. This analysis again treats the cell distribution as static and ignores the possibility that neutrophils may randomly encounter the target.

Wave initiation with Hill function activation, cells in 1D, and diffusion in 1D

For cells in 1D with diffusion in 1D, the scenario described above can be described with the following equation of motion:

(79) ct=D2cr2+aρh2Θ[ri-|r|]+aρh2Θ[|r|-ri]cncn+Cthn.

One can non-dimensionalize Equation (79) by dividing all the concentration scales by Cth, dividing all the length scales by lc=h2CthD/aρ, and dividing all the time scales by τc=h2Cth/aρ, thusly arriving at

(80) (c/Cth)(t/τc)=2(c/Cth)(r/lc)2+Θ[ri-|r|lc]+Θ[|r|-rilc](c/Cth)n(c/Cth)n+1,

which shows that lc=h2CthD/aρ is the relevant length scale and τc=h2Cth/aρ is the relevant time scale. (In the n limit of Heaviside activation, τc=D/v2 and lc=D/v.) As with Heaviside activation, we refer to the initiation time tinit as the time at which c(ri,tinit)=Cth. To find tinit, we numerically solve Equation (79) using the methods discussed in Appendix 5: Asymptotic wave dynamics with Hill function activation. This gives the relationship of tinit/τc as a function of ri/lc and n shown in Appendix 8—figure 1A. As seen in Appendix 8—figure 1A, even low-order (n=1,2,3,5,10) Hill functions exhibit relatively large (compared to τc) initiation times for rilc.

Appendix 8—figure 1
Diffusive wave initiation with Hill function activation.

(A) Wave initiation for cells in 1D with diffusion in 1D. Here, the initiation time tinit is normalized by the characteristic time scale τc=h2Cth/aρ and the initial signaling colony size ri is normalized by the characteristic length scale lc=(h2CthD/aρ)1/2. Yellow data points are individual simulations with dashed black guides to the eye connecting the yellow points. The solid yellow line corresponds to Heaviside-like activation and is reproduced from the main text. (B) Same as A, but for cells in 2D and diffusion in 2D. Here, τc=hCth/aρ and lc=(hCthD/aρ)1/2. (C) Same as A, but for cells in 3D and diffusion in 3D. Here, τc=Cth/aρ and lc=(CthD/aρ)1/2. We note the value of ri=3lc, the value below which waves fail to initiate for Heaviside-like activation.

Wave initiation with Hill function activation, cells in 2D, and diffusion in 2D

Next, we study the initiation dynamics above with cells and diffusion in two dimensions. Here, the dynamics are governed by the following equation of motion:

(81) ct=D(2cr2+1rcr)+aρhΘ[ri-r]+aρhΘ[r-ri]cncn+Cthn.

Now that we are considering the initiation dynamics, we must include terms like D(c/r)/r, which we could previously neglect in our asymptotic analysis of cells in 2D with diffusion in 2D. The curvature of the initial signaling colony matters when calculating initiation times. Note that Equation (81) can be non-dimensionalized in the same spirit as Equation (80) with characteristic length scale lc=hCthD/aρ and characteristic time scale τc=hCth/aρ. As with cells and diffusion in 1D, we numerically solve Equation (81) using the methods discussed in Appendix 5: Asymptotic wave dynamics with Hill function activation to find tinit such that c(ri,tinit)=Cth. Here again, we see that even low-order (n=2,3,5,10) Hill functions exhibit relatively large (compared to τc) initiation times for rilc.

Wave initiation with Hill function activation, cells in 3D, and diffusion in 3D

Finally, we study signaling wave initiation properties of a 3D environment by studying wave initiation with cells and diffusion in 3D. To do so, we numerically solve the 3D analog of Equation (81),

(82) ct=D(2cr2+2rcr)+aρΘ[ri-r]+aρΘ[r-ri]cncn+Cthn,

which can be non-dimensionalized in the same spirit as Equation (80), but with characteristic length scale lc=CthD/aρ and characteristic time scale τc=Cth/aρ. We solve Equation (82) using the methods discussed in Appendix 5: Asymptotic wave dynamics with Hill function activation to find tinit such that c(ri,tinit)=Cth. This gives the numerically determined relationship of tinit as a function of ri and n shown in Appendix 8—figure 1C. We see that even low-order (n=2,3) Hill functions exhibit relatively large (compared to τc) initiation times for rilc. Larger yet Hill functions (n=5,10) can lead to very large initiation times (compared to τc) even for rilc.

In fact, n>3 activation functions can result in initiation failures when aρri2/3DCth1. To see that this is the case, we treat Equation (82) in the steady state (c/t=0) and use a perturbative analysis, assuming cn/(cn+Cthn)1, in which case cn/(cn+Cthn)(c/Cth)n. In such a situation, we can write c(r)c0(r)+c1(r) as the sum of a dominant contribution c0 that is generated by cells within ri and satisfies

(83) 0=D(2c0r2+2rc0r)+aρΘ[ri-r]

and a small correction c1 that is generated by cells beyond ri and obeys

(84) 0=D(2c1r2+2rc1r)+aρc0nc0n+CthnΘ[r-ri]D(2c1r2+2rc1r)+aρ(c0Cth)nΘ[r-ri].

We can solve Equation (83) directly to arrive at:

(85) c0(r<ri)=aρ2D(ri2-r2/3),c0(r>ri)=aρri33Dr

where the form of c0(r>ri) is reminiscent of solving for the potential of a uniformly charged sphere in electrostatics (Berg and Purcell, 1977). If aρri2/3DCth=ϵ1, we can calculate the perturbation c1. By Equation (84), c1 obeys

(86) 0=D(2c1r2+2rc1r)+aρ(c0Cth)nΘ[r-ri]=D(2c1r2+2rc1r)+aρ(ϵrir)nΘ[r-ri]

so that

(87) c1(r<ri)=aρri2ϵnD(n-2),c1(r>ri)=aρϵnD(n-3)[ri3Dr+r2n-2(rir)n].

For c1 to be a sensible perturbative correction, we require it to be positive (since we are adding source terms to c0 to get c1) and much smaller than c0. This is the case when n>3. In this limit, it is smaller than c0 by roughly a factor of ϵn – a very small correction. Thus, n>3 activation functions can give steady-state concentration profiles that do not trigger waves in a three-dimensional diffusive environment.

In Appendix 8—figure 1C, we can indeed see that small n activation functions show less appreciable increases in the initiation time as ri decreases.

Appendix 9

Sensitivity of the information front to fit parameters

As an example of the Green’s function method described in Materials and methods, we have plotted various information wave fronts in Appendix 9—figure 1. These wave fronts assume different values of the diffusion constant D and the threshold concentration Cth is fit to give the experimentally observed wave initiation time in Reátegui et al., 2017. As we can see from the plot, there is a small range of values for D for which one can construct an information wave front that agrees with the data (Appendix 9—figure 1 A). These values of D are consistent with the diffusion constant of small molecules like LTB4. Values of D differing significantly from this range give information wave fronts that differ significantly from the experimentally observed wave front.

Appendix 9—figure 1
Wavefront fitting of the simple diffusion and relay models.

(A) Plot showing the information waves for various choices of diffusion constant (black) overlaid on the experimental chemotactic index data from Reátegui et al., 2017 (color plot). Here, we take ri=100 µm, although the target in the experiment is a smaller, oblong object. The size of the target has no effect on the convex shape of the information wave front. The black line is reproduced from the main text and has D=1.25×10-10 m2/s and asymptotic wave speed v1.7 µm/s (the threshold concentration is Cth/aρ=2/πv). This diffusion constant is consistent with a small molecule like LTB4, and the resulting information wave dynamics can be made to fit the information wave front – both the initiation time of 200 s and the transient dynamics. Other choices of parameters (green: D=1.8×10-10 m2/s, v2.3 µm/s, navy: D=0.8×10-10 m2/s, v1.3 µm/s) give information wave fronts that are also roughly consistent with the experimental data. (B) However, with a much larger (D=10-9 m2/s, red) or smaller (D=10-11 m2/s, dashed) diffusion constant, an information wave with the correct initiation time does not have the correct transient dynamics. The wave speeds for these larger and smaller diffusion constants are, respectively, v11 µm/s and v0.3 µm/s. (C) Simple diffusion models of various diffusion constants overlaid atop data from Reategui et al. in which the the same swarming assay as in A/B (but with a slightly smaller 80 µm target) was utilized but with neutrophils whose LTB4 receptors (BLT1/2) had been blocked. Here, we see a 250 s offset before the propagation of a slow-moving diffusive front. A simple diffusive model does not capture this offset well, but does accurately capture the concave shape of the front (which contrasts with the convex shape of the front propagated by a relay). (D) Same as in C but with a 250 s offset in the theoretical curves. These curves fit the observed chemotactic index dynamics fairly well.

Physiological relevance of these fit parameters

With the fit values of Cth/aρ=3.67×105 s/m and D=1.25×10-10 m2/s reported in the main text and the empirical cell density of ρ=(1/50 µm)2, we conclude that Cth/a1.5×1014 s/m3.

Under similar conditions with the swarming assay, (Reátegui et al., 2017) report neutrophil LTB4 production in a 3 mm thick extracellular medium to be approximately 1 pg per 100 microliters per hour. With cells at a density of ρ=(1/50 µm)2 and a 3 mm thick extracellular medium, this corresponds to a production rate of a40 molecules/second/cell. Thus, combining this production rate with the above, we yield Cth500 pM.

This value is within the range of the measured BLT1 receptor affinity for LTB4, which is reported to be approximately 0.1 – 2 nM (Yokomizo, 2015).

Appendix 10

Simple diffusion model

In the main text, we showed the qualitative differences between a signaling relay, in which cells emit one type of signaling molecule in response to the local concentration of the same molecule, and a simple diffusive signaling model, in which cells within some volume signal surrounding cells, which do not participate in the signaling at all. Here, we explicitly calculate some of the properties of a simple diffusion model. We consider cells in 2D, with a region of cells of size ri at z=0 in which the cells emit diffusible signaling molecules at rate a. In equation form,

(88) ct=D(2cr2+1rcr+2cz2)+aρδ(z)Θ[ri-r]

describes the concentration in space and time. To calculate concentrations, one can either propagate this equation directly or, as we do in the main text, integrate Green’s functions in a manner similar to that described above. For cells in 2D with diffusion in 3D (assuming a semi-infinite environment), this gives

(89) c(r,z=0,t)=aρ0tdT0ridRG2,3(r,0;R,-T)

as the concentration at z=0. One can take gradients according to c/r.

In Appendix 9—figure 1, we show that these dynamics produce information wave fronts that are different from the information wave fronts of diffusive relays. Moreover, the information wave fronts of simple diffusive theories are inconsistent (Appendix 9—figure 1C) with the observed chemotactic index dynamics catalogued by Reátegui et al., 2017 for neutrophils with blocked LTB4 (BLT1/2) receptors. With a time delay of ∼250 s, the information wave fronts of simple diffusive models can be made consistent with the observed information wave fronts. As there is great ambiguity about the signaling molecules that govern recruitment of BLT1/2-blocked neutrophils, we have shown this is true for a large range (two orders of magnitude) of diffusion constants.

In the limit rri, the signaling colony looks like a point source, meaning that Equation (89) can be simplified according to 𝑑RG2,3ri2e-r2/4DT/(4πD3/2T3/2). Thus, we get that

(90) rri:c(r,z=0,t)aρri22πD3/20tdTe-r2/4DTT3/2=aρri22Drerfc(r24Dt)1/2

which, in the limit of r2Dt, gives

(91) rri,Dt:c(r,z=0,t)aρri2Drerfc(r24Dt)1/2aρri2r2tπDe-r2/4Dt

which shows that the concentration profiles of simple diffusive models indeed have very shallow, Gaussian (with 1/r2 adjustments) tails.

Similarly, for cells in 2D with diffusion in 2D,

(92) c(r,t)=aρh0t𝑑T0ri𝑑RG2,2(r,0;R,-T)

describes the concentrations. In the same limits (rri,Dt), we get

(93) rri,Dt:c(r,t)aρri2thr2e-r2/4Dt

as the concentration. Here again, we see that the concentration profiles are shallow Gaussian with 1/r2 adjustments. We plot the resulting gradients in this thin extracellular medium limit against the gradients from a comparable relay model in Appendix 10—figure 1B.

Appendix 10—figure 1
Comparison of relay and simple diffusion for cells and diffusion in 2D.

All plots are for D=10-10 m2/s, v=2 µm/s, and ri=4D/v. (A) Information fronts for relay (black) and simple diffusion (gray) models. The information wave travels like t for simple diffusion and vt for the relay. The threshold concentration for the simple diffusion model is chosen such that the its information front intersects the relay model at t=20D/v2. Thus, hCth,relay/aρ=25 s while hCth,diff./aρ0.25 s. (B) Snapshots of radial gradients at various times for both the relay (black) and simple diffusion (gray) models. The dashed vertical lines indicate the location of the wave fronts. For short times (left) the relay’s information wave front lags behind the simple diffusion model’s information front, though it later catches up (middle) and passes it (right). At all times, for cells just inside the wave front, the relay model creates gradients that are orders of magnitude larger than does simple diffusion.

Appendix 11

Quantifying the effects of chemotaxis

To understand the effects that chemotaxis has on our model, we consider cells in 2D. The results below can be adopted to study cells in 3D or 1D, although the dimensionality of the cells has no effect on the asymptotic signaling wave speed.

We consider the same signaling motif as in the main text – that cells emit a diffusible molecule with rate a once the local concentration of the same molecule exceeds Cth – but now consider a time-varying density. Our model is a coarse-grained one; we study the case of cells moving toward the origin (radially inward) with a mean speed u once the local concentration exceeds Cth. This is a toy model of neutrophil chemotaxis and is, of course, an approximation because the mean radial speed will depend on – among many other factors – the strengths of the gradients the cells use for chemotaxis. In full, for cells in any number of dimensions and within this model,

(94a) c(𝐫,t)t=D2c+aρ(𝐫,t)Θ[c-Cth]
(94b) ρ(𝐫,t)t=(ur^ρ)Θ[c(𝐫,t)-Cth]

where r^ is the unit vector pointing radially outward. For cells in 2D,

(95a) c(r,z,t)t=D(2cr2+1rcr+2cz2)+aρ(r,t)δ(z)(Θ[c-Cth]Θ[r-ri]+Θ[ri-r])
(95b) ρ(r,t)t=ur(rρ)rΘ[c-Cth]Θ[r-ri]

which are the coupled equations we will study going forward. Note that we have included a source term for cells within the initial signaling colony of radius ri.

We again assume a signaling wave propagates at outward with speed v. In this case, the cell density beyond the target is described by:

(96) ρ(ri<r<vt,t)=ρ01+ut/r(1+u/v)2

which one can derive by assuming an outward propagating wave with speed v, a group of inward chemotaxing cells with speed u, and an initially uniform density of cells ρ0. To do so, we consider the signaling wave passing a cell at radius R and time t-T. At a later time t, the cell initially at R has moved inward a distance uT. Thus, the density at r=R-uT and t is ρ(r,t)ρ0R/r=ρ0(1+ut/r)/(1+u/v). Integrating this density and demanding conservation of cell number gives Equation (96).

Effect on asymptotic wave speed relationships

Before numerically solving Equations (95a) and (95b) to show the precise effects chemotaxis has on the concentration profiles, concentration profiles, and information wave fronts, we calculate the effect it has on the asymptotic wave speed.

To do so, we first show that only cells within D/v of the wave front contribute to the concentration at the wave front. To contribute to the concentration at the wave front, you need to be within about a diffusion length of it. If the wave front passed a time t ago, that means being within δrDt. However, we know that t=δr/v, so δrD/v – the characteristic length scale of diffusive waves.

As only cells within D/v of the wave front contribute to concentration at the wave front, we are considering cell densities on the order of:

(97) ρ(r=vt-D/v,t)=ρ01+ut/r(1+u/v)2ρ01+ut/(vt-D/v)(1+u/v)2.

In the asymptotic regime of vtD/v, we get

(98) ρρ0/(1+u/v)

meaning the density of cells that contribute to the wave front propagation is approximately constant. Therefore, we may modify the analysis that lead to Equations (24) and (29) to get two new asymptotic equations for the wave speed:

(99) hD/v:Cth=aρ0Dhv2(1+u/v)


(100) hD/v:Cth=2aρ0πv(1+u/v).

For neutrophils and the information wave front presented in the main text (reproduced in Appendix 11—figure 1), 1+u/v1+(0.3μm/s)/(2μm/s)=1.15 and the effect of chemotaxis on the asymptotic dynamics is small.

Appendix 11—figure 1
Effects of chemotaxis on wavefront location and concentration profiles.

(A) Information wave fronts for cell signaling relays with (navy) and without (black) chemotaxis. The information wave fronts are overlaid on the experimental chemotactic index data from (color plot) (Reátegui et al., 2017). The black curve is reproduced from the main text. Both models can account for the observed information wave fronts by fitting two parameters: the signaling molecule diffusivity, D and the threshold concentration, Cth. (B/C) Concentration profiles (B) and gradients (C) generated by the signaling relay models in A. The wave front is indicated in all panels by the dashed line, and the concentration profiles and gradients are plotted at the times such that the threshold concentration is equal to the concentration at the wave front. When one accounts for chemotaxis, the concentration profiles near the target steepen relative to models without chemotaxis.

Effect on transient dynamics

To study the effect of chemotaxis on the transient dynamics of neutrophil swarming, we utilize a modified version of our Green’s function method reported in Materials and methods to propagate Equations (95a) and Equations (95b). One can do so using the same algorithm described previously, but with the following modifications with rc(t) the information wave front:

  • For radii r at time t satisfying ri<r<rc(t), find the time t*(r,t) at which the neutrophils at radius r at time t began chemotaxing inward. This time satisfies the relationship rc[t*]-r=u(t-t*).

  • The density at r is therefore given by the ratio of rc[t*] to r with an additional factor of one plus the ratio of the inward advection speed to the outward-propagating wave speed at rc[t*] and t*:

    (101) ρ(r,t)=ρ0rc[t*]/r1+u(rc[t*]t)-1.

    This is analogous to Equation (96) and reduces to Equation (96) in the asymptotic limit of rc(t)=vt.

  • We assume that once the cells reach the target edge, they pack inward at a maximum density, in units of the cell diameter dc, of ρmax=1/dc2. For the experiments we discuss (Reátegui et al., 2017), this means that ρmax10ρ0.

With all these adjustments, and using the reported value of u20 µm/s in Reátegui et al., 2017, we arrive at the navy information wave front in Appendix 11—figure 1A. This curve is a fit by eye to the experimental information wave front and has fit parameters of D=1.5×10-10 m2/s and v=1.73 µm/s, corresponding to a threshold concentration of Cth/aρ0=2πv(1+u/v)3.07×105 s/m. For reference, the black curve in Appendix 11—figure 1A is the information wave front from Figure 4 of the main text, for which D=1.25×10-10 m2/s and Cth/aρ0=2πv(1+u/v)3.67×105 s/m. Thus, including chemotaxis only negligibly affects our fit values.

To compare the two models in Appendix 11—figure 1A, we plot both the concentration profile (Appendix 11—figure 1B) and the concentration gradient (Appendix 11—figure 1C) for a given critical radius (the dashed line in Appendix 11—figure 1A). When one accounts for chemotaxis, the concentration profile near the target steepens relative to a model with a stationary cell distribution. We can therefore see that chemotaxis itself can lead to steeper concentration profiles, though the model we have explored here only accounts for an average inward drift of the cells.

Data availability

The only dataset we analyze or generate is present and available in Reategui (2017). PMID:29057147. Code for the figures is available at https://github.com/pdieterle/diffWavePropAndInit (copy archived at https://archive.softwareheritage.org/swh:1:rev:f8d9feffd57d05f47c8c14c6d9850643b2858d0a/).


    1. Kessler DA
    2. Levine H
    (1993) Pattern formation in Dictyostelium via the dynamics of cooperative biological entities
    Physical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 48:4801–4804.
    1. Luther R
    Raumliche fortplanzung chimischer reaktionen
    Z Elektrochem 12:596–600.

Article and author information

Author details

  1. Paul B Dieterle

    Department of Physics, Harvard University, Cambridge, United States
    Formal analysis, Investigation, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8129-7456
  2. Jiseon Min

    Department of Molecular and Cellular Biology, Harvard University, Cambridge, United States
    Formal analysis, Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Daniel Irimia

    BioMEMS Resource Center and Center for Surgery, Innovation and Bioengineering, Department of Surgery, Massachusetts General Hospital, Boston, United States
    Conceptualization, Supervision, Methodology, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7347-2082
  4. Ariel Amir

    John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, United States
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Investigation, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2611-0139


National Science Foundation (MRSEC DMR 14-20570)

  • Ariel Amir

Kavli Foundation

  • Ariel Amir

National Science Foundation (CAREER 1752024)

  • Ariel Amir

Fannie and John Hertz Foundation (Paul M. Young Fellowship)

  • Paul B Dieterle

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


We thank Richard Murray, James Parkin, Justin Bois, Mikhail Shapiro, David Nelson, Taekjip Ha, Johanna Dickmann, and Jenny Sheng for helpful discussions; Felice Frankel for help with figure preparation; and the reviewers of this manuscript for their thorough reading and helpful comments. We acknowledge support from the NSF through MRSEC DMR 14–20570 and the Kavli Foundation. PBD is supported by the Paul M Young Fellowship through the Fannie and John Hertz Foundation. AA acknowledges support from NSF CAREER 1752024.

Version history

  1. Received: August 4, 2020
  2. Accepted: December 3, 2020
  3. Accepted Manuscript published: December 4, 2020 (version 1)
  4. Version of Record published: January 4, 2021 (version 2)


© 2020, Dieterle et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 2,812
  • 371
  • 16

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Paul B Dieterle
  2. Jiseon Min
  3. Daniel Irimia
  4. Ariel Amir
Dynamics of diffusive cell signaling relays
eLife 9:e61771.

Share this article


Further reading

    1. Developmental Biology
    2. Physics of Living Systems
    Raphaël Clément

    Geometric criteria can be used to assess whether cell intercalation is active or passive during the convergent extension of tissue.

    1. Computational and Systems Biology
    2. Physics of Living Systems
    Taegon Chung, Iksoo Chang, Sangyeol Kim
    Research Article

    Locomotion is a fundamental behavior of Caenorhabditis elegans (C. elegans). Previous works on kinetic simulations of animals helped researchers understand the physical mechanisms of locomotion and the muscle-controlling principles of neuronal circuits as an actuator part. It has yet to be understood how C. elegans utilizes the frictional forces caused by the tension of its muscles to perform sequenced locomotive behaviors. Here, we present a two-dimensional rigid body chain model for the locomotion of C. elegans by developing Newtonian equations of motion for each body segment of C. elegans. Having accounted for friction-coefficients of the surrounding environment, elastic constants of C. elegans, and its kymogram from experiments, our kinetic model (ElegansBot) reproduced various locomotion of C. elegans such as, but not limited to, forward-backward-(omega turn)-forward locomotion constituting escaping behavior and delta-turn navigation. Additionally, ElegansBot precisely quantified the forces acting on each body segment of C. elegans to allow investigation of the force distribution. This model will facilitate our understanding of the detailed mechanism of various locomotive behaviors at any given friction-coefficients of the surrounding environment. Furthermore, as the model ensures the performance of realistic behavior, it can be used to research actuator-controller interaction between muscles and neuronal circuits.