Dynamics of diffusive cell signaling relays
Abstract
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 fasttraveling 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\approx {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 quasi1D tubes (Cheng and Ferrell, 2018; Chang and Ferrell, 2013; Nolet et al., 2020), in quasi2D 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 wellunderstood 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 quasi1D 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 singlecomponent 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 wellknown scaling laws can be drastically altered. For example, cells confined to two (or one) dimensions with signaling molecule diffusion in three (or two) dimensions give rise to a diffusive wave whose speed has no dependence on D: a wave driven by diffusion whose speed does not depend on the rate of diffusion. In contrast to the dramatic effect of system dimensionality, these scaling laws are insensitive to many biological details, including the functional form of cellular activation — the dependence of signaling molecule emission rate on the local concentration. We additionally account for other phenomena – molecule decay, pulsed emission, and the discreteness of cells – that do affect the asymptotic wave dynamics; in so doing, we provide an intuitive rubric for determining under what conditions these effects alter the wave propagation speed.
In our studies of wave initiation, we systematically examine under what conditions a group of cells can trigger the formation of a diffusive wave. Here again, our results provide predictive relationships between biophysical inputs and the resulting dynamics, which are at once dramatically affected by dimensionality and largely insensitive to the details of activation and cellular uptake.
Finally, we show that neutrophil swarming experiments (Reátegui et al., 2017) display dynamics consistent with our model. In this context, our results elucidate a potential design principle of diffusive relays: they create large concentration gradients. Whereas simple diffusion of a signaling molecule from a central source creates a shallow concentration profile that falls off like $\mathrm{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 singlecomponent relays contrast with the pulse trains of chemotactic cues observed in, for example, Dictyostelium discoideum (Kessler and Levine, 1993; Pálsson and Cox, 1996).
Results
Model construction
We begin by considering a static group of cells uniformly distributed in two dimensions – for example, atop a solid surface – and described by an area density $\rho $ (Figure 1A). We assume a cell at position $\mathbf{\mathbf{r}}$ senses the local concentration of a signaling molecule, $c(\mathbf{\mathbf{r}},t)$, and participates in the emission at a concentrationdependent 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{\mathbf{r}},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 nonlinear function of c. It can include uptake, release, and cellinduced 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 twodimensional plane, one can study the signaling dynamics of cells in a onedimensional channel or a threedimensional 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(\mathbf{\mathbf{r}},t)=c(r,z,t)$ can be described by an outwardtraveling wave of the form $c(r,z,t)=c(\stackrel{~}{r}=rvt,z)$ (Fisher, 1937; Kolmogorov et al., 1937; Tanaka et al., 2017). Here, $\stackrel{~}{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(\stackrel{~}{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 \stackrel{~}{r})/r$ term due to the dominance of $v\partial c/\partial \stackrel{~}{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(\stackrel{~}{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^{2}/s and a wave speed of $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 $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 outofplane 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 twodimensional 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 longpulse, longdecay 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 nonsource 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^{2}/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 nondimensionalizing 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 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) for $\stackrel{~}{r}>0$ and $\stackrel{~}{r}<0$, then matching boundary conditions at $\stackrel{~}{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 FisherKolmogorov 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^{2}/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 timedependent dynamics of the source term at the wave front, a phenomenon that breaks our assumption that all cells obey the same timeindependent 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:
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 semiinfinite 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^{2}/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 singlecell activation.
It is worth reflecting on the fact that some system geometries give a wave whose speed is diffusion constantindependent. 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 zdimension 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 Dindependent 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 dimensiondependent wave initiation phenomena (Foerster et al., 1989; Weise and Panfilov, 2011); our task here is to study the dimensiondependent dynamics of concentration build up in singlecomponent 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 twodimensional diffusive environments, a continuously emitting source leads to a diverging concentration throughout the space. However, in threedimensional diffusive environments, a continuously emitting source gives rise to a steadystate 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 twodimensional environments; meanwhile, a colony in a threedimensional environment may fail to initiate a diffusive wave. This is indeed what we observe.
In Figure 3C, we demonstrate this dramatic dimensiondependence in the case of switchlike 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 dimensiondependent dynamics of concentration buildup dictate that a signaling wave which will always initiate in one and twodimensional environments requires a critical initial colony size in 3D.
Because the signaling wave always initiates in one and twodimensional 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 twodimensional 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 positivevalued 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 fullscale 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 Hilllike activation function $f(c)={c}^{n}/({c}^{n}+{C}_{\text{th}}^{n})$ of order $n\ge 2$, the above results for switchlike 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 switchlike activation. These results highlight the role of spatial degrees of freedom in determining the wave initiation dynamics and stability.
Application to neutrophil swarming and gradient generation
With a firm understanding of the diffusive wave and initiation dynamics, we now turn our sights to understanding a specific model system: neutrophil swarming. In beautiful work across several organisms (Lämmermann et al., 2013; Isles et al., 2019; Reátegui et al., 2017), experimentalists have observed striking behavior: an acute injury or infection can elicit rapid, highly directive motion of neutrophils – the most prevalent white blood cells – toward the site of the injury or infection. These experiments have demonstrated that a lipid small molecule called leukotriene B4 (LTB4) – along with many larger, slowerdiffusing proteins (Reátegui et al., 2017) – governs the longrange 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 antiinflammatory 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, humanderived 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 $\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 $\u27e8\mathrm{cos}\theta \u27e9$ 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}\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 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\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^{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\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, $d=50$ µm, satisfies $vd/4D\approx 0.17\ll 1$. The cell thickness $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 fieldofview and longer timecourse 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 $C}_{\text{th}}\approx 500\text{}\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 wellknown, a burstlike 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\approx 100\mu $ µm), one could probe the thin extracellular medium limit of $h\ll D/v$ with microfluidic chambers of tens of microns in height. Similarly, with mmscale 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 dimensionalitydependent 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 fasttravelling wave, then how do the cells collectively turn off response? One possibility is that signaling pathways in neutrophil swarm resolution – for instance, those involving LXA4 (Reátegui et al., 2017) production and emission – work by a similar relay mechanism; it is also possible that LTB4 production is governed by other 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; $dTdRa\rho {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 $\rho $ 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}_{\text{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\rho dTdR{G}_{n,m}({r}_{c}(t),t;R,T)$ at ${r}_{c}(t)$; the sum of all these contributions must be equal to ${C}_{\text{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 finegrained meshing over the outofplane 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 inplane sources; the Green’s functions appropriately keep track of the outofplane 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 $\mathrm{\Delta}t\ll D/{v}^{2}$ (we use $\mathrm{\Delta}t=D/10{v}^{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)$.
An explicit implementation of this method is provided at github./pdieterle/diffWavePropAndInit (Dieterle, 2020; copy archived at swh:1:rev:f8d9feffd57d05f47c8c14c6d9850643b2858d0a).
Appendix 1
Model setup
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{\mathbf{r}},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 pointlike cells homogeneously distributed in a twodimensional plane, we will write the source term as $a\rho \delta (z)\mathrm{\Theta}[c{C}_{\text{th}}]$. Here, $\rho $ is a twodimensional cell density, $\delta (.)$ is the Dirac delta function, and z is the outofplane dimension. Therefore, the cell density $\rho $ which appears in all of the subsequent discussion will be a threedimensional density if the cells are distributed in three dimensions, a twodimensional 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 outofplane dimensions is justified when the cell size is small compared to the inherent length scale of the diffusive wave, $D/v$.
Appendix 2
Asymptotic wave ansatz
We start out by seeking to understand what dynamical properties a system described by Equation (14) has at large times. To study these dynamics, we need an inspired guess for what the dynamics will look like. We imagine that when a small volume of cells starts signaling its neighbors – and those neighbors start signaling their neighbors – that a reasonable guess for the dynamics of such a signaling relay is an outward propagating wave with speed v. We define r as the distance from the center of the outward propagating wave. At long times for a uniform cell density, the shape information wave front will obey radial symmetry and our ansatz becomes $c(\mathbf{\mathbf{r}},t)=c(\mathbf{\mathbf{r}}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 $\stackrel{~}{r}=rvt$ (we call this $\stackrel{~}{x}=xvt$ for cells in one dimension) which defines the distance to the wave front. Note that $\stackrel{~}{r}<0$ means we are inside the wave front while $\stackrel{~}{r}>0$ means we are beyond it. With these definitions, $\partial c/\partial \stackrel{~}{r}=\partial c/\partial r$ and $\partial c/\partial t=v\partial c/\partial \stackrel{~}{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 outofplane 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 \stackrel{~}{r}$ dominates terms like $D\left(\partial c/\partial \stackrel{~}{r}\right)/r$ and we can ignore the latter (Tanaka et al., 2017). By construction, our ansatz says that $c(\stackrel{~}{r}=0)\equiv {C}_{\text{th}}$, meaning that $\mathrm{\Theta}[c{C}_{\text{th}}]=\mathrm{\Theta}[\stackrel{~}{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 nondimensionalizing, 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 onedimensional 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 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 nondimensional 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}\u27f9v\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 twodimensional diffusive environment (the thin extracellular medium limit) and an effectively threedimensional environment (the thick extracellular medium limit).
Cells in 2D, diffusion in 2D: the thin 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 zdependence and concentrations get normalized by h. Here,
is our asymptotic governing equation.
For both $\stackrel{~}{r}<0$ and $\stackrel{~}{r}>0$, Equation (21) reduces to two straightforwardtosolve 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 $\stackrel{~}{r}\to \mathrm{\infty}$ and that the concentration profile only blow up linearly as $\stackrel{~}{r}\to \mathrm{\infty}$. Physically, these demands are justified as follows: the concentration as $\stackrel{~}{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 $\stackrel{~}{r}\to \mathrm{\infty}$ can grow at most linearly because the cells a distance $\stackrel{~}{r}$ from the wave front have only been emitting for a time $\stackrel{~}{r}/v$. Combining these asymptotic boundary conditions with the demand that $c(\stackrel{~}{r})$ be continuous and have a continuous first derivative at $\stackrel{~}{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.
Cells in 2D, diffusion in 3D: the thick extracellular medium limit
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 semiinfinite 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}}\mathrm{exp}({z}^{2}/2{H}^{2})$ and
Next, we take a partial Fourier transform of Equation (26) with k and $C(\stackrel{~}{r},k)$ the Fourier partners of z and $c(\stackrel{~}{r},z)$, respectively. Here, we choose $C(\stackrel{~}{r},k)\equiv \frac{1}{\sqrt{2\pi}}{\int}_{\mathrm{\infty}}^{\mathrm{\infty}}{e}^{ikz}c(\stackrel{~}{r},z)$. This gives:
which is another pair of piecewise, straightforwardtosolve, 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 realspace concentrations and can be replaced with 1. This is equivalent to having chosen $\delta (z)$ to describe the outofplane 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 $\stackrel{~}{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 $\stackrel{~}{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 nondimensionalize 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).
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\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; $\stackrel{~}{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 selfconsistency relationship for ${C}_{\text{th}}$ by inverse Fourier transforming the analogs of Equations (31a) and Equations (31b). The value of ${C}_{\text{th}}=c(\stackrel{~}{x}=0,y=z=0)$ diverges as $H\to 0$.
Appendix 3
Pulsed emission and decay
In this section, we consider pulsed emission and decay of the signaling molecule. These scenarios are relevant for signaling pathways in, for example, Dictyostelium (Pálsson and Cox, 1996; Noorbakhsh et al., 2015; Kessler and Levine, 1993) and E. coli (Parkin and Murray, 2018), in which intracellular dynamics produce a pulselike 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 $\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 $\stackrel{~}{r}<v\tau $ (Region I) or $\stackrel{~}{r}>0$ (Region III) and $a{\rho}_{0}$ for $v\tau <\stackrel{~}{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}=2{C}_{\text{th}}/a\rho $) – the concentration profile for $\stackrel{~}{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
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\stackrel{~}{r}/D$ becomes comparable to ${v}^{2}\tau /D$, at which point ${e}^{{v}^{2}\tau /Dv\stackrel{~}{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 wavelike solution for shorter pulses.
We note that a timed pulsed emission considered here is formally equivalent to cells signaling until the local concentration exceeds $c(\stackrel{~}{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).
Continuous emission plus decay with cells in 2D and diffusion in 2D
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,
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 decayfree counterpart (Appendix 3—figure 1A). For $\gamma \ll {v}^{2}/D$, Equation (41) gives the decayfree relationship. And – as with pulsed emission – the wave speed goes to zero, this time when $\gamma \to 1/2{\tau}_{c}$ where ${\tau}_{c}=h{C}_{\text{th}}/a\rho $ (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 realspace concentration and the following selfconsistency relationship for the wave speed:
which simplifies to ${C}_{\text{th}}=2a\rho /\pi v$ for $\tau \gg D/{v}^{2}$.
We can compare $\tau $ to the characteristic time ${\tau}_{c}=D{(\pi {C}_{\text{th}}/2a\rho )}^{2}$, which is equal to $D/{v}^{2}$ in the limit of continuous emission. Numerical solution of Equation (45) reveals that there is no selfconsistent solution to Equation (45) until about $\tau \approx 1.53{\tau}_{c}$, at which point $v\approx 0.6\times 2a\rho /\pi {C}_{\mathrm{t}\mathrm{h}};\tau \approx 1.53{\tau}_{c}$ is larger than the minimum initiation time of ${t}_{\text{min,3D}}=4{\tau}_{c}/\pi $ (Appendix 3—figure 2 C, Section Appendix 7: Initiation dynamics). Thus, an initial pulse of length $4{\tau}_{c}/\pi <\tau <1.53{\tau}_{c}$ from cells within an initial signaling radius r_{i} can cause neighboring cells to exceed ${C}_{\text{th}}$, but cannot trigger a wavelike solution asymptotically.
Continuous emission plus decay with cells in 2D and diffusion in 3D
We can add signaling molecule decay to the embedded system dynamics by adding a term of $\gamma C(\stackrel{~}{r},k)$ to Equation (27) with $\gamma $ the signaling molecule decay rate. Going through the same exercise yields:
with
as the parameter relationship obtained after an inverse Fourier transform at $z=0$ and $x=0$. This gives a profile that propagates as a pulse (Appendix 3—figure 1B).
In the limit $\gamma \ll {v}^{2}/D$, we recover the familiar expression ${C}_{\text{th}}=2a\rho /\pi v$. Again, as when we accounted for signaling molecule decay with cells in 2D and diffusion in 2D, the wave speed approaches zero, but with $\gamma ={(\pi /4)}^{2}/{\tau}_{c}$ where ${\tau}_{c}=D{(\pi {C}_{\text{th}}/2a\rho )}^{2}$ (Appendix 3—figure 2 D).
Appendix 4
Finite extracellular medium
We consider now what happens when one does not lie in the extreme cases of an extracellular medium of thickness $h\gg D/v$ or $h\ll D/v$. To examine the case of arbitrary thickness h, we turn to the method of images (Appendix 4—figure 1A).
Our boundary conditions for the extracellular medium require the concentration to obey $\partial c/\partial z=0$ at $z=0,h$ – signaling molecules cannot diffuse through the boundaries. By invoking the uniqueness theorem, we know that if we can find an arrangement of 'image cells' – each emitting diffusible signaling molecule at rate a – that satisfies these boundary conditions, then this arrangement of cells gives the unique solution for the concentration profile inside the extracellular medium. In our case, to satisfy the boundary condition above, we have image cells at $z=\pm 2jh$ for $j=1,2,\mathrm{\dots}$.
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 burstlike emission by a single pointlike source at $(R,2jh,T)$ is given by the Green’s function:
so the Green’s function of a single pointlike source in a finitethickness extracellular medium is (Appendix 4—figure 1A):
We assume a traveling wave solution at speed v meaning that the concentration ${C}_{\text{th}}=c(r=vt,z=0,t\to \mathrm{\infty})$ is, for a density of cells $\rho $ emitting with rate a, given by:
with the substitutions $\stackrel{~}{t}=tT$ and $\stackrel{~}{R}=Rvt$. This yields:
for $x=v\stackrel{~}{R}/4D$ and $\mathrm{E}\mathrm{i}[.]$ the exponential integral function. For $h\gg D/v$, Equation (51) reduces to Equation (29) because the term ${\sum}_{j=1}^{\mathrm{\infty}}\mathrm{\cdots}\approx 0$. Meanwhile, for $h\ll 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 $a\rho h/\pi D$, we arrive at a relationship between a nondimensionalized threshold concentration, $\pi {C}_{\text{th}}D/a\rho h$, and a nondimensionalized 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\ll D/v$) and thick ($h\gg D/v$) extracellular medium limits.
Appendix 5
Asymptotic wave dynamics with Hill function activation
Numerical solutions show traveling waves
As shown in the main text, making the change from a Heaviside function source term to an ordern Hill function ($\mathrm{\Theta}[c{C}_{\text{th}}]\to {c}^{n}/({c}^{n}+{C}_{\text{th}}^{n})$) 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\ge 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.
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}_{\mathrm{\Theta}}=2$ µm/s with the threshold concentration determined by Equations (29) and (24) with $v\to {v}_{\mathrm{\Theta}}$. In this way, we could compare v_{n} – the wave speed given by the ordern Hill function – to the Heaviside wave speed ${v}_{\mathrm{\Theta}}$. For our numerics, we imposed a maximum size step of $D/10{v}_{\mathrm{\Theta}}$ for the spatial dimension where the cells live, a maximum step size of $D/30{v}_{\mathrm{\Theta}}$ for the added spatial dimension (when modeling diffusion in 2D), and a maximum time step of $D/10{v}_{\mathrm{\Theta}}^{2}$. 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}_{\text{th}}$. We simulated times ${t}_{\text{max}}\le 25D/{v}_{\mathrm{\Theta}}^{2}$ using Mathematica’s ‘NDSolveValue’ function and, when modeling diffusion in 2D, replaced the delta function in Equation (14) with a Gaussian of width $D/10{v}_{\mathrm{\Theta}}$. As noted earlier, since the width of this Gaussian is sufficiently smaller than D/v, the numerical results will be close to the delta function limit; indeed, this substitution gives a wave speed (according to an inverse Fourier transform of Equation (28b) at $\stackrel{~}{r},z=0$) of $v/{v}_{\mathrm{\Theta}}\approx 0.95$, in close correspondence with the Delta function wave speed. The initiating colony is of size ${r}_{i}=2D/{v}_{\mathrm{\Theta}}$, 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}_{\mathrm{\Theta}}^{2}\le t\le 25D/{v}_{\mathrm{\Theta}}^{2}$. As shown in Appendix 5—figure 1, for $n\ge 2$, the wave speeds are very close to ${v}_{\mathrm{\Theta}}$ with significant deviation only for $n=1$. Even the concentration profiles are in good agreement with the Heaviside solution for $n\ge 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 and Kolmogorov et al., 1937. To see this, we first consider a modified version of Equation (16c) with ordern 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 $\mathrm{\Theta}[c{C}_{\text{th}}]$ into a simple function of $\stackrel{~}{r}:\mathrm{\Theta}[c{C}_{\mathrm{t}\mathrm{h}}]=\mathrm{\Theta}[\stackrel{~}{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 nonlinear source term. We note that Equation (53) is, in the limit $n\to \mathrm{\infty}$, equivalent to a relay with Heaviside activation.
However, there is a distinct advantage to the new source term: it is onetoone in c. Thus, we may make the substitution $f=\frac{{c}^{n}}{{c}^{n}+{C}_{\text{th}}^{n}}\u27f9c={C}_{\text{th}}{\left(\frac{f}{1f}\right)}^{1/n}$, 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\to \{0,1\}$, a term $D\frac{{\partial}^{2}f}{{\partial}^{2}\stackrel{~}{r}}$, and a term $v\frac{\partial f}{\partial \stackrel{~}{r}}$. 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)=\frac{\partial f}{\partial x}$, which allows us to make the substitution $\frac{{\partial}^{2}f}{\partial {x}^{2}}=\frac{\partial F}{\partial x}=\frac{\partial f}{\partial x}\frac{\partial F}{\partial f}=F\frac{\partial F}{\partial f}$. 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 nonlinear 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\gg {C}_{\text{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}/\left({c}^{n}+{C}_{\text{th}}^{n}\right)$, in these regions.
Next, we turn to the traditional method of examining the $f\to 0$ limit. As $F\to 0$, Equation (56) becomes, to lowest order in f and assuming $F\approx \lambda {f}^{\beta}$,
where $\lambda >0$. One can only obtain a selfconsistency relationship between v and ${C}_{\text{th}}/a\rho $ if $\beta =21/n$. Otherwise, ${f}^{21/n\beta}$ diverges or goes to zero. With this choice, Equation (57) becomes:
where in Equation (58b) we have taken the $f\to 0$ limit. To get a wave speed, v, out of Equation (58a), we demand that the quantity under the square root be nonnegative, which ensures that $\lambda >0$ is a real number as assumed. This means $v\ge 2\sqrt{a\rho D/{C}_{\text{th}}}=2{v}_{\mathrm{\Theta}}$ – 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\sqrt{a\rho D/{C}_{\text{th}}}$ is what one finds after numerically solving the 1D dynamics (Appendix 5—figure 1A). In contrast, for $n>1$, this method yields no such wave speed bound.
Appendix 6
Assessing the validity of a continuum analysis
Next, we consider the validity of a continuum analysis like Equation (16b) for studying asymptotic wave dynamics. To do so, we compare our continuum wave speed relationships, Equations (24) and (29), to a simple model of discrete cells on a lattice in 1D (Appendix 6—figure 1A). We refer to this as the discrete lattice model, which has been studied previously and is the subject of ongoing work (Keener, 2000; Dieterle P and Amir A, 2020. Manuscript in preparation). We briefly discuss the results of this lattice model below, and show that it agrees with the continuous model when the separation between cells, d, is much less than the characteristic length $D/v$. This heuristic also holds for cells in two and three dimensions, and for cells scattered randomly according to a Poisson process.
To start, we first calculate the threedimensional 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}_{\bullet ,m}(x,t)$. For diffusion in $m=2$ dimensions, we will consider a semiinfinite environment in order to recapitulate Equation (29), which holds for cells in one (two) dimensions with diffusion in a semiinfinite twodimensional (threedimensional) space. These relationships are:
with $\text{erfc}[.]$ the complementary error function and $\text{Ei}[.]$ the exponential integral. We have assumed onedimensional diffusion takes place in a narrow channel with crosssectional 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}_{\text{th}}$, they participate in the signaling molecule emission at a rate a. If the resulting wave speed is v, then a cell at a distance $\stackrel{~}{x}=jd$ from the wave front has been emitting for a time $jd/v$. That cell then creates a concentration ${c}_{j,m}={c}_{\bullet ,m}(jd,jd/v)$ at $(\stackrel{~}{x},t)=(0,0)$. The full concentration at the wave front is ${C}_{\text{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 selfconsistency relationships between the threshold concentration, ${C}_{\text{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\ll 1$, the sums in these relationships are wellapproximated by an integral over j from $j=0$ to $\mathrm{\infty}$. In this limit, with $\rho =1/d$, Equation (60a) becomes ${C}_{\text{th}}=a\rho D/{h}^{2}{v}^{2}$, the onedimensional analog of Equation (24); similarly, Equation (60b) simplifies to ${C}_{\text{th}}=2a\rho /\pi hv$, the onedimensional analog of Equation (29). (One can turn Equation (60b) into an integral from $j=0$ to $\mathrm{\infty}$. 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\ll 4D/v$. We have shown the approach to the continuous theory limit in Appendix 6—figure 1B/C.
Appendix 7
Initiation dynamics
In this section, we demonstrate the initiation time relationships discussed in the main text using Green’s function integration. To do so, we write down the Green’s functions ${G}_{n,m}(r,t;R,T)$ describing diffusion of molecules in m dimensions released by cells in n dimensions at $(R,T)$ and measured by cells at $(r,t)$. For $n\ne m$, we assume a semiinfinite environment. For cells in 1D and diffusion in 1D, we assume a narrow channel of width h in both dimensions perpendicular to the channel. For cells in 1D or 2D and diffusion in 2D, we assume an extracellular medium of thickness h. We calculate the Green’s functions for cells in two dimensions by integrating over a ring of diffusive sources at radius R; we calculate the Green’s functions for cells in three dimensions by integrating over a shell of diffusive sources at radius R. Below, ${I}_{0}[.]$ is the zeroth IBessel function and sinh $[.]$ is the hyperbolic sine function. The Green’s functions are
Onebyone, we study the initiation time for these systems by studying the selfconsistency relationship for the threshold concentration ${C}_{\text{th}}$ and initiation time ${t}_{\text{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}_{\text{th}}$. At ${t}_{\text{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}_{\text{init}}$ and initial signaling colony size r_{i} for various system dimensionalities when ${r}_{i}\gg D/v$ or ${r}_{i}\ll D/v$. The full relationships of ${t}_{\text{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}\gg D{t}_{i}$, we get that ${t}_{i}=2D/{v}^{2}$ by using the asymptotic relationship Equation (24). This is the minimum initiation time, ${t}_{\text{min},1,1}$.
and it tells us that ${r}_{i}\gg D{t}_{i}$ is equivalent to ${r}_{i}\gg 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}\gg D/v$ initiation time.
In the opposite limit – ${r}_{i}\ll D/v$ – we can Taylor Expand Equation (63) and get
thus validating our equations in the main text.
Initiation with cells in 1D and diffusion in 2D
Next, we consider the selfconsistency equation Equation (62) with $n=1,m=2$. As before, we first consider the limit of ${r}_{i}\gg Dv$ and recover (through Equation (29)):
while for ${r}_{i}\ll Dv$,
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}\gg D/v$, such an analysis gives a minimum initiation time of
In the opposite limit of ${r}_{i}\ll 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}_{\text{init}}\to \mathrm{\infty}$ which gives us a maximum concentration ${C}_{\text{max},2,3}$ at $r={r}_{i}$ of:
Thus, by Equation (29), if ${r}_{i}<D/v$, ${C}_{\text{max},2,3}<{C}_{\text{th}}$ and the signaling wave cannot initiate. Examining Equation (62) for ${r}_{i}\approx D/v$ reveals
while
Initiation with cells in 3D and diffusion in 3D
Finally, we consider initiation with cells in 3D. Again, we consider the limit ${t}_{\text{init}}\to \mathrm{\infty}$ in Equation (62) to get:
Thus, waves do not initiate below a critical r_{i}. However, here, the critical value is ${r}_{i}=\sqrt{3}D/v$. As with 1D cells/diffusion and 2D cells/diffusion, we recover
in the limit ${r}_{i}\gg D/v$.
Wave initiation with discrete cells
In the main text, we claimed that the qualitative wave initiation findings hold even in systems with discrete cells. Here, we show this by considering the extreme case of a single, pointlike cell at the center of an initiating colony of size r_{i}. We consider cells and diffusion in the same number of dimensions, m. In an mdimensional diffusive environment, the concentration created by this source – which emits at a rate a – at a radius r_{i} and time t will be given by:
For $m=3$, this integral is bounded from above by $a/4\pi D{r}_{i}$ as $t\to \mathrm{\infty}$ while for $m=1,2$, the concentration diverges as $t\to \mathrm{\infty}$. Thus we can see that, as in the continuum theory of the main text, the wave will always initiate for $m=1,2$ but has a critical radius for $m=3$.
However, this is not the only qualitative similarity between the discrete and continuum cases. For $m=1,2$, we obtain
which we can manipulate to a more familiar form by considering the ${r}_{i}\ll D/v$ and $t\gg D/{v}^{2}$ limit, in which case we yield:
Making the substitution of $c({r}_{i},t)={C}_{\text{th}}$, setting $t={t}_{\text{init}}$, multiplying both sides of the equations by ${r}_{i}^{m}$, and writing ${C}_{\text{th}}=a\rho D/{v}^{2}\sim aD/{r}_{i}^{m}$, we arrive at:
exactly the results of Equations (65) and (69).
Thus, not only the rough phenomenology (of $m=3$ being distinct from $m=1,2$ in that exhibits a critical radius), but also the specific scalings of initiation time hold when accounting for the discreteness of cells.
Appendix 8
Wave initiation with Hill function activation
We now characterize the effect of Hill functionlike activation on the signaling wave initiation, as we have done already for asymptotic dynamics. To do so, we consider a simple situation: cells within a volume of size r_{i} signal with some rate a while neighboring cells outside of the initial signaling volume (i.e., with $r>{r}_{i}$) signal with a concentrationdependent rate of $a{c}^{n}/({c}^{n}+{C}_{\text{th}}^{n})$. These calculations give us an idea of how sensitive the initiation dynamics are to the details of cell activation. As we will show, the initiation dynamics with Hill function activation are a good approximation of the initiation dynamics for Heaviside activation when for relatively small n. One can imagine that such a situation may be relevant in, for example, the neutrophil swarming experiments presented in the main text (Reátegui et al., 2017). Here, cells in direct contact with a foreign protein begin signaling their neighbors, which respond to the presence of the signaling molecule by participating in the emission themselves. This analysis again treats the cell distribution as static and ignores the possibility that neutrophils may randomly encounter the target.
Wave initiation with Hill function activation, cells in 1D, and diffusion in 1D
For cells in 1D with diffusion in 1D, the scenario described above can be described with the following equation of motion:
One can nondimensionalize Equation (79) by dividing all the concentration scales by ${C}_{\text{th}}$, dividing all the length scales by ${l}_{c}=\sqrt{{h}^{2}{C}_{\text{th}}D/a\rho}$, and dividing all the time scales by ${\tau}_{c}={h}^{2}{C}_{\text{th}}/a\rho $, thusly arriving at
which shows that ${l}_{c}=\sqrt{{h}^{2}{C}_{\text{th}}D/a\rho}$ is the relevant length scale and ${\tau}_{c}={h}^{2}{C}_{\text{th}}/a\rho $ is the relevant time scale. (In the $n\to \mathrm{\infty}$ limit of Heaviside activation, ${\tau}_{c}=D/{v}^{2}$ and ${l}_{c}=D/v$.) As with Heaviside activation, we refer to the initiation time ${t}_{\text{init}}$ as the time at which $c({r}_{i},{t}_{\text{init}})={C}_{\text{th}}$. To find ${t}_{\text{init}}$, we numerically solve Equation (79) using the methods discussed in Appendix 5: Asymptotic wave dynamics with Hill function activation. This gives the relationship of ${t}_{\text{init}}/{\tau}_{c}$ as a function of ${r}_{i}/{l}_{c}$ and n shown in Appendix 8—figure 1A. As seen in Appendix 8—figure 1A, even loworder ($n=1,2,3,5,10$) Hill functions exhibit relatively large (compared to ${\tau}_{c}$) initiation times for ${r}_{i}\ll {l}_{c}$.
Wave initiation with Hill function activation, cells in 2D, and diffusion in 2D
Next, we study the initiation dynamics above with cells and diffusion in two dimensions. Here, the dynamics are governed by the following equation of motion:
Now that we are considering the initiation dynamics, we must include terms like $D(\partial c/\partial r)/r$, which we could previously neglect in our asymptotic analysis of cells in 2D with diffusion in 2D. The curvature of the initial signaling colony matters when calculating initiation times. Note that Equation (81) can be nondimensionalized in the same spirit as Equation (80) with characteristic length scale ${l}_{c}=\sqrt{h{C}_{\text{th}}D/a\rho}$ and characteristic time scale ${\tau}_{c}=h{C}_{\text{th}}/a\rho $. As with cells and diffusion in 1D, we numerically solve Equation (81) using the methods discussed in Appendix 5: Asymptotic wave dynamics with Hill function activation to find ${t}_{\text{init}}$ such that $c({r}_{i},{t}_{\text{init}})={C}_{\text{th}}$. Here again, we see that even loworder ($n=2,3,5,10$) Hill functions exhibit relatively large (compared to ${\tau}_{c}$) initiation times for ${r}_{i}\ll {l}_{c}$.
Wave initiation with Hill function activation, cells in 3D, and diffusion in 3D
Finally, we study signaling wave initiation properties of a 3D environment by studying wave initiation with cells and diffusion in 3D. To do so, we numerically solve the 3D analog of Equation (81),
which can be nondimensionalized in the same spirit as Equation (80), but with characteristic length scale ${l}_{c}=\sqrt{{C}_{\text{th}}D/a\rho}$ and characteristic time scale ${\tau}_{c}={C}_{\text{th}}/a\rho $. We solve Equation (82) using the methods discussed in Appendix 5: Asymptotic wave dynamics with Hill function activation to find ${t}_{\text{init}}$ such that $c({r}_{i},{t}_{\text{init}})={C}_{\text{th}}$. This gives the numerically determined relationship of ${t}_{\text{init}}$ as a function of r_{i} and n shown in Appendix 8—figure 1C. We see that even loworder ($n=2,3$) Hill functions exhibit relatively large (compared to ${\tau}_{c}$) initiation times for ${r}_{i}\ll {l}_{c}$. Larger yet Hill functions ($n=5,10$) can lead to very large initiation times (compared to ${\tau}_{c}$) even for ${r}_{i}\approx {l}_{c}$.
In fact, $n>3$ activation functions can result in initiation failures when $a\rho {r}_{i}^{2}/3D{C}_{\text{th}}\ll 1$. To see that this is the case, we treat Equation (82) in the steady state ($\partial c/\partial t=0$) and use a perturbative analysis, assuming ${c}^{n}/({c}^{n}+{C}_{\text{th}}^{n})\ll 1$, in which case ${c}^{n}/({c}^{n}+{C}_{\text{th}}^{n})\approx {\left(c/{C}_{\text{th}}\right)}^{n}$. In such a situation, we can write $c(r)\approx {c}_{0}(r)+{c}_{1}(r)$ as the sum of a dominant contribution c_{0} that is generated by cells within r_{i} and satisfies
and a small correction c_{1} that is generated by cells beyond r_{i} and obeys
We can solve Equation (83) directly to arrive at:
where the form of ${c}_{0}(r>{r}_{i})$ is reminiscent of solving for the potential of a uniformly charged sphere in electrostatics (Berg and Purcell, 1977). If $a\rho {r}_{i}^{2}/3D{C}_{\text{th}}=\u03f5\ll 1$, we can calculate the perturbation c_{1}. By Equation (84), c_{1} obeys
so that
For c_{1} to be a sensible perturbative correction, we require it to be positive (since we are adding source terms to c_{0} to get c_{1}) and much smaller than c_{0}. This is the case when $n>3$. In this limit, it is smaller than c_{0} by roughly a factor of ${\u03f5}^{n}$ – a very small correction. Thus, $n>3$ activation functions can give steadystate concentration profiles that do not trigger waves in a threedimensional diffusive environment.
In Appendix 8—figure 1C, we can indeed see that small n activation functions show less appreciable increases in the initiation time as r_{i} decreases.
Appendix 9
Sensitivity of the information front to fit parameters
As an example of the Green’s function method described in Materials and methods, we have plotted various information wave fronts in Appendix 9—figure 1. These wave fronts assume different values of the diffusion constant D and the threshold concentration ${C}_{\text{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.
Physiological relevance of these fit parameters
With the fit values of ${C}_{\text{th}}/a\rho =3.67\times {10}^{5}$ s/m and $D=1.25\times {10}^{10}$ m^{2}/s reported in the main text and the empirical cell density of $\rho =(1/50$ µm)^{2}, we conclude that ${C}_{\text{th}}/a\approx 1.5\times {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 $\rho =(1/50$ µm)^{2} and a 3 mm thick extracellular medium, this corresponds to a production rate of $a\approx 40$ molecules/second/cell. Thus, combining this production rate with the above, we yield ${C}_{\text{th}}\approx 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).
Appendix 10
Simple diffusion model
In the main text, we showed the qualitative differences between a signaling relay, in which cells emit one type of signaling molecule in response to the local concentration of the same molecule, and a simple diffusive signaling model, in which cells within some volume signal surrounding cells, which do not participate in the signaling at all. Here, we explicitly calculate some of the properties of a simple diffusion model. We consider cells in 2D, with a region of cells of size r_{i} at $z=0$ in which the cells emit diffusible signaling molecules at rate a. In equation form,
describes the concentration in space and time. To calculate concentrations, one can either propagate this equation directly or, as we do in the main text, integrate Green’s functions in a manner similar to that described above. For cells in 2D with diffusion in 3D (assuming a semiinfinite environment), this gives
as the concentration at $z=0$. One can take gradients according to $\partial c/\partial r$.
In Appendix 9—figure 1, we show that these dynamics produce information wave fronts that are different from the information wave fronts of diffusive relays. Moreover, the information wave fronts of simple diffusive theories are inconsistent (Appendix 9—figure 1C) with the observed chemotactic index dynamics catalogued by Reátegui et al., 2017 for neutrophils with blocked LTB4 (BLT1/2) receptors. With a time delay of ∼250 s, the information wave fronts of simple diffusive models can be made consistent with the observed information wave fronts. As there is great ambiguity about the signaling molecules that govern recruitment of BLT1/2blocked neutrophils, we have shown this is true for a large range (two orders of magnitude) of diffusion constants.
In the limit $r\gg {r}_{i}$, the signaling colony looks like a point source, meaning that Equation (89) can be simplified according to $\int \mathit{d}R{G}_{2,3}\approx {r}_{i}^{2}{e}^{{r}^{2}/4DT}/\left(4\sqrt{\pi}{D}^{3/2}{T}^{3/2}\right)$. Thus, we get that
which, in the limit of ${r}^{2}\gg Dt$, gives
which shows that the concentration profiles of simple diffusive models indeed have very shallow, Gaussian (with $1/{r}^{2}$ adjustments) tails.
Similarly, for cells in 2D with diffusion in 2D,
describes the concentrations. In the same limits ($r\gg {r}_{i},\sqrt{Dt}$), we get
as the concentration. Here again, we see that the concentration profiles are shallow Gaussian with $1/{r}^{2}$ adjustments. We plot the resulting gradients in this thin extracellular medium limit against the gradients from a comparable relay model in Appendix 10—figure 1B.
Appendix 11
Quantifying the effects of chemotaxis
To understand the effects that chemotaxis has on our model, we consider cells in 2D. The results below can be adopted to study cells in 3D or 1D, although the dimensionality of the cells has no effect on the asymptotic signaling wave speed.
We consider the same signaling motif as in the main text – that cells emit a diffusible molecule with rate a once the local concentration of the same molecule exceeds ${C}_{\text{th}}$ – but now consider a timevarying density. Our model is a coarsegrained one; we study the case of cells moving toward the origin (radially inward) with a mean speed u once the local concentration exceeds ${C}_{\text{th}}$. This is a toy model of neutrophil chemotaxis and is, of course, an approximation because the mean radial speed will depend on – among many other factors – the strengths of the gradients the cells use for chemotaxis. In full, for cells in any number of dimensions and within this model,
where $\widehat{r}$ is the unit vector pointing radially outward. For cells in 2D,
which are the coupled equations we will study going forward. Note that we have included a source term for cells within the initial signaling colony of radius r_{i}.
We again assume a signaling wave propagates at outward with speed v. In this case, the cell density beyond the target is described by:
which one can derive by assuming an outward propagating wave with speed v, a group of inward chemotaxing cells with speed u, and an initially uniform density of cells ${\rho}_{0}$. To do so, we consider the signaling wave passing a cell at radius R and time $tT$. At a later time t, the cell initially at R has moved inward a distance $uT$. Thus, the density at $r=RuT$ and t is $\rho (r,t)\sim {\rho}_{0}R/r={\rho}_{0}(1+ut/r)/(1+u/v)$. Integrating this density and demanding conservation of cell number gives Equation (96).
Effect on asymptotic wave speed relationships
Before numerically solving Equations (95a) and (95b) to show the precise effects chemotaxis has on the concentration profiles, concentration profiles, and information wave fronts, we calculate the effect it has on the asymptotic wave speed.
To do so, we first show that only cells within $\approx D/v$ of the wave front contribute to the concentration at the wave front. To contribute to the concentration at the wave front, you need to be within about a diffusion length of it. If the wave front passed a time t ago, that means being within $\delta r\approx \sqrt{Dt}$. However, we know that $t=\delta r/v$, so $\delta r\approx D/v$ – the characteristic length scale of diffusive waves.
As only cells within $\approx D/v$ of the wave front contribute to concentration at the wave front, we are considering cell densities on the order of:
In the asymptotic regime of $vt\gg D/v$, we get
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
For neutrophils and the information wave front presented in the main text (reproduced in Appendix 11—figure 1), $1+u/v\approx 1+(0.3\mu \text{m/s})/(2\mu \text{m/s})=1.15$ and the effect of chemotaxis on the asymptotic dynamics is small.
Effect on transient dynamics
To study the effect of chemotaxis on the transient dynamics of neutrophil swarming, we utilize a modified version of our Green’s function method reported in Materials and methods to propagate Equations (95a) and Equations (95b). One can do so using the same algorithm described previously, but with the following modifications with ${r}_{c}(t)$ the information wave front:
For radii r at time t satisfying ${r}_{i}<r<{r}_{c}(t)$, find the time ${t}^{*}(r,t)$ at which the neutrophils at radius r at time t began chemotaxing inward. This time satisfies the relationship ${r}_{c}\left[{t}^{*}\right]r=u(t{t}^{*})$.
The density at r is therefore given by the ratio of ${r}_{c}\left[{t}^{*}\right]$ to r with an additional factor of one plus the ratio of the inward advection speed to the outwardpropagating wave speed at ${r}_{c}[{t}^{*}]$ and ${t}^{*}$:
(101) $\rho (r,t)={\rho}_{0}\frac{{r}_{c}\left[{t}^{*}\right]/r}{1+u{\left(\frac{\partial {r}_{c}\left[{t}^{*}\right]}{\partial t}\right)}^{1}}.$This is analogous to Equation (96) and reduces to Equation (96) in the asymptotic limit of ${r}_{c}(t)=vt$.
We assume that once the cells reach the target edge, they pack inward at a maximum density, in units of the cell diameter d_{c}, of ${\rho}_{\text{max}}=1/{d}_{c}^{2}$. For the experiments we discuss (Reátegui et al., 2017), this means that ${\rho}_{\text{max}}\approx 10{\rho}_{0}$.
With all these adjustments, and using the reported value of $u\approx 20$ µm/s in Reátegui et al., 2017, we arrive at the navy information wave front in Appendix 11—figure 1A. This curve is a fit by eye to the experimental information wave front and has fit parameters of $D=1.5\times {10}^{10}$ m^{2}/s and $v=1.73$ µm/s, corresponding to a threshold concentration of ${C}_{\text{th}}/a{\rho}_{0}=\frac{2}{\pi v(1+u/v)}\approx 3.07\times {10}^{5}$ s/m. For reference, the black curve in Appendix 11—figure 1A is the information wave front from Figure 4 of the main text, for which $D=1.25\times {10}^{10}$ m^{2}/s and ${C}_{\text{th}}/a{\rho}_{0}=\frac{2}{\pi v(1+u/v)}\approx 3.67\times {10}^{5}$ s/m. Thus, including chemotaxis only negligibly affects our fit values.
To compare the two models in Appendix 11—figure 1A, we plot both the concentration profile (Appendix 11—figure 1B) and the concentration gradient (Appendix 11—figure 1C) for a given critical radius (the dashed line in Appendix 11—figure 1A). When one accounts for chemotaxis, the concentration profile near the target steepens relative to a model with a stationary cell distribution. We can therefore see that chemotaxis itself can lead to steeper concentration profiles, though the model we have explored here only accounts for an average inward drift of the cells.
Data availability
The only dataset we analyze or generate is present and available in Reategui (2017). PMID:29057147. Code for the figures is available at https://github.com/pdieterle/diffWavePropAndInit (copy archived at https://archive.softwareheritage.org/swh:1:rev:f8d9feffd57d05f47c8c14c6d9850643b2858d0a/).
References

LTB4 is a signalrelay molecule during neutrophil chemotaxisDevelopmental Cell 22:1079–1091.https://doi.org/10.1016/j.devcel.2012.02.003

Spatial waves of advance with bistable dynamics: cytoplasmic and genetic analogues of allee effectsThe American Naturalist 178:E48–E75.https://doi.org/10.1086/661246

Physics of chemoreceptionBiophysical Journal 20:193–219.https://doi.org/10.1016/S00063495(77)855446

The pharmacology and pathophysiology of leukotriene B4British Medical Bulletin 39:249–254.https://doi.org/10.1093/oxfordjournals.bmb.a071828

Engineering microbial consortia: a new frontier in synthetic biologyTrends in Biotechnology 26:483–489.https://doi.org/10.1016/j.tibtech.2008.05.004

Chemical waves in cell and developmental biologyJournal of Cell Biology 217:1193–1204.https://doi.org/10.1083/jcb.201701158

The wave of advance of advantageous genesAnnals of Eugenics 7:355–369.https://doi.org/10.1111/j.14691809.1937.tb02153.x

Spatial trigger waves: positive feedback gets you a long wayMolecular Biology of the Cell 25:3486–3493.https://doi.org/10.1091/mbc.e14081306

Propagation and its failure in coupled systems of discrete excitable cellsSIAM Journal on Applied Mathematics 47:556–572.https://doi.org/10.1137/0147038

Propagation of waves in an excitable medium with discrete release sitesSIAM Journal on Applied Mathematics 61:317–334.https://doi.org/10.1137/S0036139999350810

Pattern formation in Dictyostelium via the dynamics of cooperative biological entitiesPhysical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 48:4801–4804.https://doi.org/10.1103/physreve.48.4801

Neutrophil swarming: an essential process of the neutrophil tissue responseImmunological Reviews 273:76–93.https://doi.org/10.1111/imr.12458

Study of the diffusion equation coupled to an increase in mass and its application to a problem in biologyBull Univ Mosc Ser Int Sec A A1:1–26.https://doi.org/10.1201/97803678105047

BookA Kinetic View of Statistical PhysicsCambridge University Press.https://doi.org/10.1017/CBO9780511780516

The multiple faces of leukocyte interstitial migrationSeminars in Immunopathology 36:227–251.https://doi.org/10.1007/s0028101404188

Signal propagation and failure in discrete autocrine relaysPhysical Review Letters 93:118101.https://doi.org/10.1103/PhysRevLett.93.118101

Modeling oscillations and spiral waves in Dictyostelium populationsPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 91:062711.https://doi.org/10.1103/PhysRevE.91.062711

Microscale arrays for the profiling of start and stop signals coordinating humanneutrophil swarmingNature Biomedical Engineering 1:0094.https://doi.org/10.1038/s415510170094

Epithelial cell guidance by selfgenerated EGF gradientsIntegrative Biology 4:259–269.https://doi.org/10.1039/c2ib00106c

Two distinct leukotriene B4 receptors, BLT1 and BLT2Journal of Biochemistry 157:65–71.https://doi.org/10.1093/jb/mvu078
Article and author information
Author details
Funding
National Science Foundation (MRSEC DMR 1420570)
 Ariel Amir
Kavli Foundation
 Ariel Amir
National Science Foundation (CAREER 1752024)
 Ariel Amir
Fannie and John Hertz Foundation (Paul M. Young Fellowship)
 Paul B Dieterle
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Richard Murray, James Parkin, Justin Bois, Mikhail Shapiro, David Nelson, Taekjip Ha, Johanna Dickmann, and Jenny Sheng for helpful discussions; Felice Frankel for help with figure preparation; and the reviewers of this manuscript for their thorough reading and helpful comments. We acknowledge support from the NSF through MRSEC DMR 14–20570 and the Kavli Foundation. PBD is supported by the Paul M Young Fellowship through the Fannie and John Hertz Foundation. AA acknowledges support from NSF CAREER 1752024.
Version history
 Received: August 4, 2020
 Accepted: December 3, 2020
 Accepted Manuscript published: December 4, 2020 (version 1)
 Version of Record published: January 4, 2021 (version 2)
Copyright
© 2020, Dieterle et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,812
 views

 371
 downloads

 16
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

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

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