Dynamics of diffusive cell signaling relays

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.


Introduction
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 (D » 10 À10 m 2 /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 cellsthat 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ðÀr 2 =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 dimensionsfor example, atop a solid surface -and described by an area density ( Figure 1A). We assume a cell at position r senses the local concentration of a signaling molecule, cðr; 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ðr; tÞ: where the Dirac delta function dð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 th . In such a case, the activation function f ðcÞ is well described by a Heaviside step function Q½c À C th and the concentration dynamics obey qc qt ¼ Dr 2 c þ adðzÞQ½c À C th : 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 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ðr; 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 beyondand 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 th cðr ¼ 0; z ¼ 0Þ and r ) D=v, we take Equation (2)  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, C th . 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 outwardpropagating waves. Here, D ¼ 10 À10 m 2 /s, v ¼ 2 mm/s, and hC th =a ¼ D=v 2 . Numerical simulations assume that a cell colony of size r i ¼ 4D=v (dashed vertical line) centered at the origin starts signaling t ¼ 0. (C) Numerical wave speed as measured at t ¼ 100D=v 2 (markers) agrees well with theory (Equation (5), black line) as we independently vary D (circles) and hC th =a (diamonds) relative to the panel B values (red circle and diamond).
Since we consider r ) D=v, we may ignore the Dðqc=qrÞ=r term due to the dominance of vqc=qr. 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 vdepend on the biophysical system parameters like the cell density, ; the concentration threshold, C 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 » 10 À10 m 2 /s and a wave speed of v » 1 mm/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=v » 100 mm and D=v 2 » 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 dðzÞ is valid when the cell size H ( D=v; similarly, decay of the signaling molecule can be neglected for a decay rate g ( ðD=v 2 Þ À1 while pulsed emission gives rise to the same asymptotic wave speed if the width of the pulse t satisfies t ) 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 ( 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 ( D=v, diffusion of the signaling molecule is effectively two-dimensional as we can take q 2 c=qz 2 ! 0 and dðzÞ ! 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=h, all concentrations in the problem, including C th , 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, C th , a, , and h combine to give us a single model parameter to describe the threshold concentration, hC th =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 m 2 /s). Thus, the only combination of these two parameters that will give a speed (measured in m/s) is ðaD=hC th Þ 1=2 . By this simple dimensional analysis argument, the wave speed v can only be v ¼ aðaD=hC th Þ 1=2 for some constant a. 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 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 cellinduced 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) forr>0 andr<0, then matching boundary conditions atr ¼ 0. This analysis indeed reveals that 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 ) 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 hC th =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 v~ffi ffiffi ffi D p , 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 2 /s) and hC th =a (measured in s), we add a modulation time t (measured in seconds) and decay rate g (measured in 1/s). As D is the only parameter involving a length scale, it must be that v~ffi ffiffi ffi D p even in these more complex scenarios.
By way of contrast, Vergassola et al., 2018 have shown that an unconventional scaling of v~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 ) 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: 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): C th =a (measured in s/m) and D (m 2 /s). The only combination of these parameters that will give a wave speed (measured in m/s) is a=C th . It therefore must be the case that v ¼ aa=C th with a 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 th and a single emission rate a. (For Hill function activation, a » 2=p for n ! 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. 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 m 2 /s and v ¼ 2 mm/s with C th =a ¼ 2=pv. The initial signaling colony is of size r i ¼ 4D=v (dashed vertical line). (C) Numerical wave speed as measured at t ¼ 100D=v 2 (markers) agrees well with theory (Equation (8), black line) as we independently vary D (circles) and C th =a (diamonds) relative to the panel B values (blue circle and diamond). As predicted, v is indeed D-independent in this system.
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 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 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).
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.  Figure 3. Wave initiation dynamics. (A) Schematic demonstrating wave initiation. Cells within some initial signaling volume of radius r i begin signaling at some rate a. The signaling wave is initiated when the concentration at nearby cells exceeds the threshold concentration, C th . (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 þ ðC th =cÞ n À1 . C: Initiation times for Heaviside activation, in which f ðcÞ ¼ Q½c À C th . Numerics (thick colored lines) and approximate asymptotic theory (Equations (10) and (11), thin black lines) of the initiation time's dependence on r i 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 vr i =D ( 1 and vr i =D ) 1. Similarly, for cells and diffusion in 2D, Equation (11) governs the large and small vr i =D limits. In both of these cases, the wave always initiates, but the initiation time Figure 3 continued on next page In Figure 3C, we demonstrate this dramatic dimension-dependence in the case of switch-like activation, for which f ðcÞ ¼ Q½c À C th . To find the initiation time t 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 th ¼ cðr i ; t 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 ( D=v) and large (r i ) D=v) initiating colonies is: When r i ( 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 th . Correspondingly, the initiation time increases like 1=r 2 i for small r i . Meanwhile, for r i ) D=v, the size of the initiating colony becomes irrelevant and reaches a minimum value of t min;1D , determined entirely by the characteristic time scale D=v 2 . The full dependence of t 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 ( D=v, the initiation time scales harshly as . Here again, the asymptotic theory Equation (12) is in good agreement with numerics. 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, t init ; vri D ( 1; and the initiation time for large initial signaling colony size, t init ; vri D ) 1. 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/m 2 for cells in 2D, and 1/m 3 for cells in 3D. When the diffusive and cell dimensions do not match, the environment is assumed to be semi-infinite.
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 vr 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 ¼ ffiffi ffi 3 p D=v below which the wave will not initiate. Around the critical colony size, t init diverges as ðvr i =DÞ 6 ½ðvr i = ffiffi ffi 3 p DÞ 2 À 1 À2 . If r i ) D=v, then t init again plateaus at a constant value t 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 buildup 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 oneand 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 n th Þ of order n ! 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.

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 ( (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).
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 hcos i 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 r i » 100 mm, 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 switchlike 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 u » 0:3 mm/s ( v (see Appendix 11: Quantifying the effects of chemotaxis). Thus, as mentioned above, Equation (2) effectively has two parameters: C th =a and D. Fitting these two parameters to the observed information wave front gives C th =a » 3:67 Â 10 5 s/m and D » 1:25 Â 10 À10 m 2 /s, the latter of which is consistent with the diffusion constant of a small molecule like LTB4. This implies a wave speed of v » 1:7 mm/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 mm, satisfies vd=4D » 0:17 ( 1. The cell thickness H » 10 mm indeed satisfies H ( D=v, meaning the use of the delta function to describe the cell distribution is valid. Finally, as LTB4 has a lifetime 1=g of many minutes (Bray, 1983) and D=v 2 » 40 s ( 1=g, 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 th =a ¼ 3:67 Â 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 » 40 molecules per second per cell (see Appendix 9: Sensitivity of the information front to fit parameters for details). Using the cell density of ¼ 1=d 2 ¼ ð50 mÞ À2 , we find that C th » 500 pM. This value is within the range of the  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 m 2 /s. The threshold concentrations for the relay and simple diffusion models are C th =a » 3:66 Â 10 5 s/m and 2:91 Â 10 4 s/m, respectively. The parameters for the relay Figure 4 continued on next page 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~ffi ffiffiffiffi Dt p 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 relaycalculated 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.

Discussion
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=v » 100 mm), one could probe the thin extracellular medium limit of h ( D=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 (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.
resolution -for instance, those involving LXA4 (Reátegui et al., 2017) production and emissionwork by a similar relay mechanism; it is also possible that LTB4 production is governed by other fastdiffusing 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, G n;m ðr; t; R; TÞ. These equations are enumerated in Appendix 7: Initiation dynamics; dT dR a G n;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 r c ðtÞ such that C th ¼ cðr c ðtÞ; tÞ. Thus, with an initial signaling colony of size r i , one must solve the problem: This constraint equation considers every radius at time T and, if it is less than r c ðTÞ, adds a concentration contribution of a dT dR G n;m ðr c ðtÞ; t; R; TÞ at r c ðtÞ; the sum of all these contributions must be equal to C th . 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 r i . 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 r c ðtÞ at discrete times, incrementing in steps of Dt ( D=v 2 (we use Dt ¼ D=10v 2 in the main text and Appendices, which gives convergence of the information wave front). Linear interpolation between these points defines a continuous curve r c ðtÞ.
is supported by the Paul M Young Fellowship through the Fannie and John Hertz Foundation. AA acknowledges support from NSF CAREER 1752024. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

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, C 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ðr; tÞ. In general, then, we have with Q½: 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 aQ½c À C th , 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 dðzÞQ½c À C th . Here, is a two-dimensional cell density, dð:Þ 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 onedimensional 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.

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 neighborsand 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ðr; tÞ ¼ cðr À vtrÞ withr the unit vector pointing from the origin to the wave front.
With this guess, we can define a new coordinater ¼ r À vt (we call thisx ¼ x À vt for cells in one dimension) which defines the distance to the wave front. Note thatr<0 means we are inside the wave front whiler>0 means we are beyond it. With these definitions, qc=qr ¼ qc=qr and qc=qt ¼ Àv qc=qr. For cells in 1D, we consider the y and z to be dimensions perpendicular to the line of cells with the density described by dðyÞdð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 dð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 These equations can be simplified once more by noting that we are considering asymptoticthat is, large r -dynamics. Thus, as long as v ) D=r, we can say that vqc=qr dominates terms like D qc=qr ð Þ=r and we can ignore the latter (Tanaka et al., 2017). By construction, our ansatz says that cðr ¼ 0Þ C th , meaning that Q½c À C th ¼ Q½À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 cells in 3D : 0 ¼ q 2 ðcv 2 =aDÞ qðvr=DÞ 2 þ qðcv 2 =aDÞ qðvr=DÞ þ Q½À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=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 ( 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 ( D=v (or an extracellular medium of thickness h ( 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 ( D=v, we can simplify Equation (16a) because q 2 c=qy 2 ¼ q 2 c=qz 2 ¼ 0 and dðyÞdðzÞ ! 1=h 2 . The resulting equation is the exact same as for cells in 3D but with a source term proportional to a=h 2 : For cells in 2D with an extracellular medium of thickness h ( D=v, diffusion is effectively twodimension 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 aD=v 2 , meaning that there must be a relationship along the lines of C th~a D=v 2 v~aD=C th ð Þ 1=2exactly 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 ) 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).

Cells in 2D, diffusion in 2D: the thin extracellular medium limit
For an extracellular medium of thickness h ( 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 bothr<0 andr>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 ! 0 asr ! ¥ and that the concentration profile only blow up linearly asr ! À¥. Physically, these demands are justified as follows: the concentration asr ! ¥ has to go to zero because there are no cells emitting in that region and it is far from the wave front; the concentration asr ! À¥ 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 atr ¼ 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.
Cells in 2D, diffusion in 3D: the thick extracellular medium limit Next, we consider Equation (16b) in the limit that the extracellular medium h ) 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 dðzÞ, we consider the cells to be of a thickness H such that 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Þ 1 ffiffiffiffi 2p p R ¥ À¥ e ikz cð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 ) v=D. Thus, if Hv=D ( 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 dðzÞ to describe the out-of-plane cell density.
Proceeding with e ÀH 2 k 2 =2 ! 1, one can take the inverse partial Fourier transform and arrive at C th ¼ 2a=pv: (29) Similarly, in the limit jrj ) D=v, e Àvr=D e Àvz 2 =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 timer=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=v, one can show that Equation (25) from which we can tell that the concentration scales of the problem, including Cth, are proportional to a=v, and thus that v~a=C th , 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 d ( D=v and diffusion in an environment of size h ) 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;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 ! k 2 y þ k 2 z . One can then find the concentration profiles and self-consistency relationship for C th by inverse Fourier transforming the analogs of Equations (31a) and Equations (31b). The value of C th ¼ cðx ¼ 0; y ¼ z ¼ 0Þ diverges as H ! 0.

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 t emitted once a cell exceeds the threshold concentration C th . In the moving frame, this pulse has length vt -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 whenr< À vt (Region I) orr>0 (Region III) and a 0 for Àvt <r<0 (Region II). We thus recover the following: Applying the same boundary conditions as with continuous emission, we arrive at RegionII :cðÀvt <r<0Þ ¼ aD=hv 2 À aD hv 2 e Àv 2 t =D e Àvr=D À ar=hv (37b) RegionIII :cðr>0Þ ¼ aD hv 2 ð1 À e Àv 2 t =D Þe Àvr=D (37c) which tells us a few things of interest. First -as seen in the right panel of Appendix 3-figure 1 A (here, t ¼ 2D=v 2 ¼ 2C th =a) -the concentration profile forr< À vt 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 For t ) D=v 2 , we recover the usual relationship of C th ¼ aD=hv 2 . In region 2, the profile will grow linearly as before until Àvr=D becomes comparable to v 2 t =D, at which point e Àv 2 t =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, ðaD=hC th Þ 1=2 , we have plotted v=ðaD=hC th Þ 1=2 as a function of 1=t in Appendix 3- figure 2A. We have normalized t by a characteristic time t c ¼ hC th =a, which is equal to D=v 2 for continuous emission. When t <t c , the wave speed goes to zero. There is no wave-like solution for shorter pulses.
1/4 1/2 1/2 1 0 1/2 1 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 t . We normalize v by the t ! ¥ wave speed of ðaD=C th Þ 1=2 . At t ¼ t c ¼ hC th =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 g. For g ¼ t min;2D ¼ 1=2t c , with t 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 t for square pulse emission of a signaling molecule by cells in 2D with a 3D diffusive environment. At t » 1:53t c ¼ 1:53Dð2C th =paÞ 2 (vertical dashed line), there is a minimum wave speed of v » 0:6 Â ð2a=pC th Þ (horizontal dashed line), unlike with 2D diffusion. Importantly, t » 1:53t c is longer than the minimum initiation time of 4t c =p. Thus, for values of t 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 g for cells in 2D with diffusion in 3D. At g ¼ ðp=4Þ 2 =t c , with t 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 ¼ Àvt Þ ¼ at =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 g by adding a term of Àgc to Equation (16b). In the thin extracellular medium limit, This is another piecewise set of linear differential equations, which we can solve as without decay to yield: as the concentration profiles and as our wave speed relationship. The concentration profile is flatter than its decay-free counterpart (Appendix 3- figure 1A). For g ( v 2 =D, Equation (41) gives the decay-free relationship. And -as with pulsed emission -the wave speed goes to zero, this time when g ! 1=2t c where t c ¼ hC th =a (Appendix 3- figure 2B).

Pulsed emission with cells in 2D and diffusion in 2D
Here, which we can solve to yield: with B i ðkÞ chosen such that C and its first derivative are continuous: Inverse Fourier transforming at z ¼ 0 gives the real-space concentration and the following selfconsistency relationship for the wave speed: which simplifies to C th ¼ 2a=pv for t ) D=v 2 .  ) showing non-dimensionalized wave speed (vh=2D) versus non-dimensionalized threshold concentration (pC th D=ah). In the limit of vh=2D ) 1, we recover the familiar 3D scaling law of Equation (29) (blue line). In the limit of vh=2D ( 1, we get the 2D scaling law of Equation (24) (red line).
Our boundary conditions for the extracellular medium require the concentration to obey qc=qz ¼ 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 ¼ AE2jh 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: Gðr; z ¼ 0; t; R; 2jh; TÞ ¼ e À ðrÀRÞ 2 þð2jhÞ 2 4DðtÀTÞ 2pDðt À TÞ ; so the Green's function of a single point-like source in a finite-thickness extracellular medium is (Appendix 4- figure 1A): We assume a traveling wave solution at speed v meaning that the concentration C th ¼ cðr ¼ vt; z ¼ 0; t ! ¥Þ is, for a density of cells emitting with rate a, given by: with the substitutionst ¼ t À T andR ¼ R À vt. This yields: for x ¼ vR=4D and Ei½: the exponential integral function. For h ) D=v, Equation (51) reduces to Equation (29) because the term P ¥ j¼1 Á Á Á » 0. Meanwhile, for h ( D=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 ah=pD, we arrive at a relationship between a non-dimensionalized threshold concentration, pC th D=ah, and a non-dimensionalized wave speed, vh=2D: We plot this relationship in Appendix 4-figure 1B and see that Equation (52) is an interpolation between the thin (h ( D=v) and thick (h ) D=v) extracellular medium limits.
As shown in the main text, making the change from a Heaviside function source term to an order-n Hill function (Q½c À C th ! c n =ðc n þ C n th Þ) 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 n ! 1 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 n ! 2, we see good agreement between the Heaviside theory (black lines) and the Hill function numerics (red dots). Snapshots are shown at t=t c ¼ 5; 10; 15; 20; 25 where t c ¼ h 2 C th =a equals D=v 2 Q , the characteristic time scale for Heaviside activation. Note that the x-axis is scaled by the characteristic length l c ¼ hðDC th =aÞ 1=2 , which is the length scale for Heaviside activation with a delta function source, D=v Q . In the insets, we display the wave speed for the order-n Hill function, v n , compared to the Heaviside activation wave speed, v Q ¼ ðaD=h 2 C th Þ 1=2 . Our fit to the n ¼ 1 data gives v n¼1 » 1:93v Q , but we have shown in Appendix 5: Asymptotic wave dynamics with Hill function activation that v n¼1 ¼ 2v Q , 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 n ! 2, 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 ¼ l c =10) Gaussian source. Here, the characteristic length l c ¼ phDC th =2a equals D=v Q , the length scale for Heaviside activation. The x-axis is scaled by the same quantity. Snapshots are shown at t=t c ¼ 5; 10; 15; 20; 25 with the characteristic time t c ¼ DðphC th =2aÞ 2 equal to D=v 2 Q , the time scale for Heaviside activation. In the insets, we display the wave speed for the order-n Hill function, v n , compared to the Heaviside activation wave speed, v Q ¼ 2a=pC th .
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 m 2 /s and v Q ¼ 2 mm/s with the threshold concentration determined by Equations (29) and (24) with v ! v Q . In this way, we could compare v n -the wave speed given by the order-n Hill function -to the Heaviside wave speed v Q . For our numerics, we imposed a maximum size step of D=10v Q for the spatial dimension where the cells live, a maximum step size of D=30v Q for the added spatial dimension (when modeling diffusion in 2D), and a maximum time step of D=10v 2 Q . These step sizes give convergence of the information wave fronts, which we define as the curves r c ðtÞ such that cðr c ðtÞ; tÞ ¼ C th . We simulated times t max 25D=v 2 Q 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 Q . 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) atr; z ¼ 0) of v=v Q » 0:95, in close correspondence with the Delta function wave speed. The initiating colony is of size r i ¼ 2D=v Q , 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, r c ðtÞ, and fit a line to the region of 24D=v 2 Q t 25D=v 2 Q . As shown in Appendix 5-figure 1, for n ! 2, the wave speeds are very close to v Q with significant deviation only for n ¼ 1. Even the concentration profiles are in good agreement with the Heaviside solution for n ! 2.

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 andKolmogorov et al., 1937. To see this, we first consider a modified version of Equation (16c) with order-n Hill function activation: 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 Q½c À C th into a simple function ofr : Q½c À C th ¼ Q½À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 , then plug this into Equation (53) to yield (after some rearrangement): Equation (54) looks a lot like the traditional Fisher equation, in that it has a source term that goes to zero at f ! f0; 1g, a term D q 2 f q 2r , and a term v qf qr . 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 Þ ¼ qf qx , which allows us to make the substitution q 2 f qx 2 ¼ qF qx ¼ qf qx qF qf ¼ F qF qf . 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: 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 c ) C th and that cells well beyond the wave front are not emitting at all. Thus, there is no spatial dependence on the cellular activation, f ¼ c n = c n þ C n th À Á , in these regions. Next, we turn to the traditional method of examining the f ! 0 limit. As F ! 0, Equation (56) becomes, to lowest order in f and assuming F » À lf b , where l>0. One can only obtain a self-consistency relationship between v and C th =a if b ¼ 2 À 1=n. Otherwise, f 2À1=nÀb diverges or goes to zero. With this choice, Equation (57) becomes: where in Equation (58b) we have taken the f ! 0 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 l>0 is a real number as assumed. This means v ! 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi aD=C th p ¼ 2v Q -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 ¼ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi aD=C th p 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.  To start, we first calculate the three-dimensional concentration (SI units of 1/m 3 ) 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: Ei À x 2 4Dt (59b) 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 h 2 and twodimensional 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 C th , they participate in the signaling molecule emission at a rate a. If the resulting wave speed is v, then a cell at a distancex ¼ Àjd from the wave front has been emitting for a time jd=v. That cell then creates a concentration c j;m ¼ c. ;m ðjd; jd=vÞ at ðx; tÞ ¼ ð0; 0Þ. The full concentration at the wave front is C th 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, C th ; wave speed, v; diffusion constant, D; and cell separation d: Equations (60a) and (60b) provide relationships analogous to Equations (24) and (29). In fact, in the limit vd=4D ( 1, 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 C th ¼ aD=h 2 v 2 , the one-dimensional analog of Equation (24); similarly, Equation (60b) simplifies to C th ¼ 2a=phv, 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 d ( 4D=v. We have shown the approach to the continuous theory limit in Appendix 6-figure 1B/C. G 3;3 ðr; t; R; TÞ ¼ R sinh rR=2Dðt À TÞ ½ e Àðr 2 þR 2 Þ=4DðtÀTÞ =r pDð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 C th and initiation time t init for a given initial signaling colony of size r i : where the logic here is that the signaling wave initiates when the threshold concentration at the edge of the initial signaling colony exceeds C th . At t init , 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 t init and initial signaling colony size r i for various system dimensionalities when r i ) D=v or r i ( D=v. The full relationships of t init versus r i 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 Àr i to r i ) and get a closed form relationship: In the limit where r i ) Dt i , we get that t i ¼ 2D=v 2 by using the asymptotic relationship Equation (24). This is the minimum initiation time, t min;1;1 .
and it tells us that r i ) Dt i is equivalent to r i ) D=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 r i ) D=v initiation time.
In the opposite limitr i ( D=v -we can Taylor Expand Equation (63) and get r i ( D=v : t init » ðpD=4v 2 ÞðD=vr i Þ 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 r i ) Dv and recover (through Equation (29)

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 r i ) D=v, such an analysis gives a minimum initiation time of t min;2;2 ¼ 2D=v 2 : In the opposite limit of r i ( D=v, 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 t init ! ¥ which gives us a maximum concentration C max;2;3 at r ¼ r i of: Thus, by Equation (29), if r i <D=v, C max;2;3 <C th and the signaling wave cannot initiate. Examining Equation (62) for r i » D=v reveals r i » D=v : t init » pD=16v 2 À Á ðvr i =DÞ 4 ðvr i =D À 1Þ À2 while t min;2;3 ¼ 4D=pv 2 :

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 C th 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. Here, we take r i ¼ 100 mm, 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 m 2 /s and asymptotic wave speed v » 1:7 mm/s (the threshold concentration is C th =a ¼ 2=pv). 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 m 2 /s, v » 2:3 mm/s, navy: D ¼ 0:8 Â 10 À10 m 2 /s, v » 1:3 mm/s) give information wave fronts that are also roughly consistent with the experimental data. (B) However, with a much larger (D ¼ 10 À9 m 2 /s, red) or smaller (D ¼ 10 À11 m 2 /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, v » 11 mm/s and v » 0:3 mm/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 mm 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 C th =a ¼ 3:67 Â 10 5 s/m and D ¼ 1:25 Â 10 À10 m 2 /s reported in the main text and the empirical cell density of ¼ ð1=50 mm) 2 , we conclude that C th =a » 1:5 Â 10 14 s/m 3 . 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 mm) 2 and a 3 mm thick extracellular medium, this corresponds to a production rate of a » 40 molecules/second/cell. Thus, combining this production rate with the above, we yield C th » 500 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).
ðr ¼ vt À D=v; tÞ ¼ 0 1 þ ut=r ð1 þ u=vÞ 2 » 0 1 þ ut=ðvt À D=vÞ ð1 þ u=vÞ 2 : In the asymptotic regime of vt ) D=v, we get » 0 =ð1 þ u=vÞ (98) 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: and h ) D=v : C th ¼ 2a 0 pvð1 þ u=vÞ : For neutrophils and the information wave front presented in the main text (reproduced in Appendix 11-figure 1), 1 þ u=v » 1 þ ð0:3 m=sÞ=ð2 m=sÞ ¼ 1:15 and the effect of chemotaxis on the asymptotic dynamics is small. 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, C th . (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.