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 $\rho$ (Figure 1A). We assume a cell at position $\mathbf{𝐫}$ senses the local concentration of a signaling molecule, $c(\mathbf{𝐫},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(\mathbf{𝐫},t)$:

where the Dirac delta function $\delta (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, $C_{\text{th}}$. In such a case, the activation function $f(c)$ is well described by a Heaviside step function $\mathrm{\Theta }[c-C_{\text{th}}]$ and the concentration dynamics obey

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.

Figure 1

Download asset Open asset

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(\mathbf{𝐫},t)=c(r,z,t)$ can be described by an outward-traveling wave of the form $c(r,z,t)=c(\tilde{r}=r-vt,z)$ (Fisher, 1937; Kolmogorov et al., 1937; Tanaka et al., 2017). Here, $\tilde{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 $C_{\text{th}}\equiv c(\tilde{r}=0,z=0)$ and $r\gg D/v$, we take Equation (2) and arrive at the following equation governing asymptotic behavior:

Since we consider $r\gg D/v$, we may ignore the $D(\partial c/\partial \tilde{r})/r$ term due to the dominance of $v\partial c/\partial \tilde{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(\tilde{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, $\rho$; the concentration threshold, $C_{\text{th}}$; 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/v^2$. For a small diffusing molecule with $D\approx {10}^{-10}$ m²/s and a wave speed of ${\displaystyle v\approx 1}$ µ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 ${\displaystyle D/v\approx 100}$ µm and $D/v^2\approx 100$ 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 $\delta (z)$ is valid when the cell size $H\ll D/v$; similarly, decay of the signaling molecule can be neglected for a decay rate $\gamma \ll {(D/v^2)}^{-1}$ while pulsed emission gives rise to the same asymptotic wave speed if the width of the pulse $\tau$ satisfies $\tau \gg D/v^2$. 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/4D\ll 1$. A thorough, mathematical discussion of all the above, including a demonstration of why $D/v$ and $D/v^2$ are the appropriate scales, is presented in Appendix 2: Asymptotic wave ansatz.

When the extracellular medium thickness $h\ll D/v$, diffusion of the signaling molecule is effectively two-dimensional as we can take ${\partial }^{2}c/\partial z^2\to 0$ and $\delta (z)\to 1/h$. In this limit, Equation (3) becomes

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\rho /h$, all concentrations in the problem, including $C_{\text{th}}$, are proportional to $a\rho /h$. As the non-source terms in Equations (4) and (1) are linear, the only role $a\rho /h$ serves is to set the concentration scale of the dynamics. Thus, $C_{\text{th}}$, a, $\rho$, and h combine to give us a single model parameter to describe the threshold concentration, $h{C}_{\text{th}}/a\rho$, 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 m²/s). Thus, the only combination of these two parameters that will give a speed (measured in m/s) is ${(a\rho D/h{C}_{\text{th}})}^{1/2}$. By this simple dimensional analysis argument, the wave speed v can only be $v=\alpha {(a\rho D/h{C}_{\text{th}})}^{1/2}$ for some constant $\alpha$. 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 $C_{\text{th}}$ 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 $\tilde{r}>0$ and $\tilde{r}<0$, then matching boundary conditions at $\tilde{r}=0$. This analysis indeed reveals that

while

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 $r\gg D/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 $h{C}_{\text{th}}/a\rho$ 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 $v\sim \sqrt{D}$, 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 m²/s) and $h{C}_{\text{th}}/a\rho$ (measured in s), we add a modulation time $\tau$ (measured in seconds) and decay rate $\gamma$ (measured in 1/s). As D is the only parameter involving a length scale, it must be that $v\sim \sqrt{D}$ even in these more complex scenarios.

By way of contrast, Vergassola et al., 2018 have shown that an unconventional scaling of $v\sim D^{3/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 $h\gg D/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:

Figure 2

Download asset Open asset

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\rho$. Thus, we have two independent parameters in Equation (7): $C_{\text{th}}/a\rho$ (measured in s/m) and D (m²/s). The only combination of these parameters that will give a wave speed (measured in m/s) is $a\rho /C_{\text{th}}$. It therefore must be the case that $v=\alpha a\rho /C_{\text{th}}$ with $\alpha$ 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 $C_{\text{th}}$ and a single emission rate a. (For Hill function activation, $\alpha \approx 2/\pi$ for $n\ge 2$, 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

and

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.

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 r_i 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).

Figure 3

Download asset Open asset

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 r_i (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)=\mathrm{\Theta }[c-C_{\text{th}}]$. To find the initiation time $t_{\text{init}}$, 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 r_i is equal to the threshold – $C_{\text{th}}=c(r_i,t_{\text{init}})$ – 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 r_i and the characteristic time and length scales – $D/v^2$ 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.

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

When $r_i\ll D/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 $C_{\text{th}}$. Correspondingly, the initiation time increases like $1/r_i^2$ for small r_i. Meanwhile, for $r_i\gg D/v$, the size of the initiating colony becomes irrelevant and reaches a minimum value of $t_{\text{min,1D}}$, determined entirely by the characteristic time scale $D/v^2$. The full dependence of $t_{\text{init}}$ on r_i 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 $r_i\ll D/v$, the initiation time scales harshly as

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 $v{r}_i/D$.

Lastly, we consider cells in 3D with 3D diffusion and find that there is a critical initial signaling colony size of $r_i=\sqrt{3}D/v$ below which the wave will not initiate. Around the critical colony size, $t_{\text{init}}$ diverges as ${(v{r}_i/D)}^6{[{(v{r}_i/\sqrt{3}D)}^2-1]}^{-2}$. If $r_i\gg D/v$, then $t_{\text{init}}$ again plateaus at a constant value $t_{\text{min,3D}}$ that only depends on the characteristic time scale $D/v^2$:

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 $r_i=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/v^2$ (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)=c^n/(c^n+C_{\text{th}}^n)$ of order $n\ge 2$, 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.

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 r_i) 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).

Figure 4

Download asset Open asset

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 $\theta$ 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 \theta ⟩$ 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 ${\displaystyle r_i\approx 100}$ µ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 ${\displaystyle u\approx 0.3}$ µm/s $\ll v$ (see Appendix 11: Quantifying the effects of chemotaxis). Thus, as mentioned above, Equation (2) effectively has two parameters: $C_{\text{th}}/a\rho$ and D. Fitting these two parameters to the observed information wave front gives $C_{\text{th}}/a\rho \approx 3.67\times {10}^5$ s/m and $D\approx 1.25\times {10}^{-10}$ m²/s, the latter of which is consistent with the diffusion constant of a small molecule like LTB4. This implies a wave speed of ${\displaystyle v\approx 1.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 $\gg D/v$ and the mean distance between neutrophils, ${\displaystyle d=50}$ µm, satisfies $vd/4D\approx 0.17\ll 1$. The cell thickness ${\displaystyle H\approx 10}$ µm indeed satisfies $H\ll D/v$, meaning the use of the delta function to describe the cell distribution is valid. Finally, as LTB4 has a lifetime $1/\gamma$ of many minutes (Bray, 1983) and $D/v^2\approx 40$ s $\ll 1/\gamma$, 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 $C_{\text{th}}/a\rho =3.67\times {10}^5$ 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 $a\approx 40$ molecules per second per cell (see Appendix 9: Sensitivity of the information front to fit parameters for details). Using the cell density of $\rho =1/d^2={(50\mu \text{m})}^{-2}$, we find that ${\displaystyle C_{\text{th}}\approx 500\mathrm{p}\mathrm{M}}$. 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 $\sim \sqrt{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.

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, $\rho$. 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, $C_{\text{th}}$. 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(\mathbf{𝐫},t)$. In general, then, we have

with $\mathrm{\Theta }[.]$ 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\rho \mathrm{\Theta }[c-C_{\text{th}}]$, with $\rho$ 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\rho \delta (z)\mathrm{\Theta }[c-C_{\text{th}}]$. Here, $\rho$ is a two-dimensional cell density, $\delta (.)$ is the Dirac delta function, and z is the out-of-plane dimension. Therefore, the cell density $\rho$ 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$.

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(\mathbf{𝐫},t)=c(\mathbf{𝐫}-vt\widehat{r})$ with $\widehat{r}$ the unit vector pointing from the origin to the wave front.

With this guess, we can define a new coordinate $\tilde{r}=r-vt$ (we call this $\tilde{x}=x-vt$ for cells in one dimension) which defines the distance to the wave front. Note that $\tilde{r}<0$ means we are inside the wave front while $\tilde{r}>0$ means we are beyond it. With these definitions, $\partial c/\partial \tilde{r}=\partial c/\partial r$ and $\partial c/\partial t=-v\partial c/\partial \tilde{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 $\rho \delta (y)\delta (z)$ with $\rho$ 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 $\rho \delta (z)$ with $\rho$ measured in cells per unit area; for cells in 3D, $\rho$ is measured in cells per unit volume. Assuming azimuthal symmetry in 2D and radial symmetry in 3D, we arrive at

These equations can be simplified once more by noting that we are considering asymptotic – that is, large r – dynamics. Thus, as long as $v\gg D/r$, we can say that $v\partial c/\partial \tilde{r}$ dominates terms like $D\left(\partial c/\partial \tilde{r}\right)/r$ and we can ignore the latter (Tanaka et al., 2017). By construction, our ansatz says that $c(\tilde{r}=0)\equiv C_{\text{th}}$, meaning that $\mathrm{\Theta }[c-C_{\text{th}}]=\mathrm{\Theta }[-\tilde{r}]$. This gives simplified equations according to:

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

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/v^2$ 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 $d\ll D/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 $h\ll D/v$ (or an extracellular medium of thickness $h\ll D/v$), then diffusion is effectively one (or two) dimensional. For instance, for cells confined in a very narrow one-dimensional channel of width $h\ll D/v$, we can simplify Equation (16a) because ${\partial }^{2}c/\partial y^2={\partial }^{2}c/\partial z^2=0$ and $\delta (y)\delta (z)\to 1/h^2$. The resulting equation is the exact same as for cells in 3D but with a source term proportional to $a\rho /h^2$:

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

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\rho D/v^2$, meaning that there must be a relationship along the lines of $C_{\text{th}}\sim a\rho D/v^2⟹v\sim {\left(a\rho D/C_{\text{th}}\right)}^{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 $h\gg D/v$. In this case, diffusion effectively takes place in three dimensions and

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).

For an extracellular medium of thickness $h\ll D/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,

is our asymptotic governing equation.

For both $\tilde{r}<0$ and $\tilde{r}>0$, Equation (21) reduces to two straightforward-to-solve linear ODEs. With b_i as constants that we will determine momentarily,

We can solve for the b_i by applying boundary conditions. First, we demand $c\to 0$ as $\tilde{r}\to \mathrm{\infty }$ and that the concentration profile only blow up linearly as $\tilde{r}\to -\mathrm{\infty }$. Physically, these demands are justified as follows: the concentration as $\tilde{r}\to \mathrm{\infty }$ 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 $\tilde{r}\to -\mathrm{\infty }$ can grow at most linearly because the cells a distance $-\tilde{r}$ from the wave front have only been emitting for a time $-\tilde{r}/v$. Combining these asymptotic boundary conditions with the demand that $c(\tilde{r})$ be continuous and have a continuous first derivative at $\tilde{r}=0$ allows us to stitch together the solutions in Equations (22a) and Equations (22b) to yield:

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

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.

Next, we consider Equation (16b) in the limit that the extracellular medium $h\gg D/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:

Instead of working directly with $\delta (z)$, we consider the cells to be of a thickness H such that $2\delta (z)\to \frac{1}{H}\sqrt{\frac{2}{\pi }}\exp (-z^2/2H^2)$ and

Next, we take a partial Fourier transform of Equation (26) with k and $C(\tilde{r},k)$ the Fourier partners of z and $c(\tilde{r},z)$, respectively. Here, we choose $C(\tilde{r},k)\equiv \frac{1}{\sqrt{2\pi }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{ikz}c(\tilde{r},z)$. This gives:

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

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 $k\gg v/D$. Thus, if $Hv/D\ll 1$, the term $e^{-H^2{k}^2/2}$ is irrelevant for calculating the real-space concentrations and can be replaced with 1. This is equivalent to having chosen $\delta (z)$ to describe the out-of-plane cell density.

Proceeding with $e^{-H^2{k}^2/2}\to 1$, one can take the inverse partial Fourier transform and arrive at

Similarly, in the limit $|\tilde{r}|\gg D/v$,

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 $\tilde{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$):

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\rho /v$, one can show that Equation (25) is equivalent to

from which we can tell that the concentration scales of the problem, including $Cth$, are proportional to $a\rho /v$, and thus that $v\sim a\rho /C_{\text{th}}$, exactly as required by Equation (29).

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 $d\ll D/v$ and diffusion in an environment of size $h\gg D/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,

with v the wave speed; $\tilde{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

which reduces to Equation (27) with $k^2\to k_y^2+k_z^2$. One can then find the concentration profiles and self-consistency relationship for $C_{\text{th}}$ by inverse Fourier transforming the analogs of Equations (31a) and Equations (31b). The value of $C_{\text{th}}=c(\tilde{x}=0,y=z=0)$ diverges as $H\to 0$.

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.

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

which is nothing more than three piecewise linear equations, which we stitch together as before. The source term is zero when $\tilde{r}<-v\tau$ (Region I) or $\tilde{r}>0$ (Region III) and $a{\rho }_0$ for $-v\tau <\tilde{r}<0$ (Region II). We thus recover the following:

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

which tells us a few things of interest. First – as seen in the right panel of Appendix 3—figure 1 A (here, $\tau =2D/v^2=2C_{\text{th}}/a\rho$) – the concentration profile for $\tilde{r}<-v\tau$ 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

Appendix 3—figure 1

Download asset Open asset

For $\tau \gg D/v^2$, we recover the usual relationship of $C_{\text{th}}=a\rho D/h{v}^2$. In region 2, the profile will grow linearly as before until $-v\tilde{r}/D$ becomes comparable to $v^2\tau /D$, at which point $e^{-v^2\tau /D-v\tilde{r}/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\rho D/h{C}_{\text{th}})}^{1/2}$, we have plotted $v/{(a\rho D/h{C}_{\text{th}})}^{1/2}$ as a function of $1/\tau$ in Appendix 3—figure 2A. We have normalized $\tau$ by a characteristic time ${\tau }_c=h{C}_{\text{th}}/a\rho$, which is equal to $D/v^2$ for continuous emission. When $\tau <{\tau }_c$, the wave speed goes to zero. There is no wave-like solution for shorter pulses.

Appendix 3—figure 2

Download asset Open asset

We note that a timed pulsed emission considered here is formally equivalent to cells signaling until the local concentration exceeds $c(\tilde{r}=-v\tau )=a\rho \tau /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).

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