# Dynamics of diffusive cell signaling relays

1. Ariel Amir 1. Department of Physics, Harvard University, United States
2. Department of Molecular and Cellular Biology, Harvard University, United States
3. BioMEMS Resource Center and Center for Surgery, Innovation and Bioengineering, Department of Surgery, Massachusetts General Hospital, United States
4. John A. Paulson School of Engineering and Applied Sciences, Harvard University, United States
Research Article
13 figures, 1 table and 1 additional file

## Figures

Figure 1 Asymptotic relay dynamics of cells in 2D with diffusion in 2D. (A) Schematic illustrating the diffusive relay motif. Cells (pink with purple nucleus) release a signaling molecule that diffuses (blue clouds). They do so when the local concentration exceeds a threshold, Cth. This gives rise to a diffusive wave with wave speed v. (B) Snapshot concentration profiles. Asymptotic theory (Equation (6), black lines) and numerical simulation of Equation (2) (red dots, details of the numerical methods can be found in Materials and methods) are in good agreement and show outward-propagating waves. Here, D=10-10 m2/s, v=2 µm/s, and h⁢Cth/a⁢ρ=D/v2. Numerical simulations assume that a cell colony of size ri=4⁢D/v (dashed vertical line) centered at the origin starts signaling t=0. (C) Numerical wave speed as measured at t=100⁢D/v2 (markers) agrees well with theory (Equation (5), black line) as we independently vary D (circles) and h⁢Cth/a⁢ρ (diamonds) relative to the panel B values (red circle and diamond).
Figure 2 Asymptotic relay dynamics with cells in 2D and diffusion in 3D. (A) Schematic of cells (pink with purple nucleus) performing a diffusive relay in which signaling molecules (blue clouds) can diffuse out-of-plane. Here, such relays give rise to a diffusion-constant-independent wave speed, v. (B) Snapshot concentration profiles of the signaling molecule show good agreement between numerical simulation of Equation (2) (blue dots, details of numerical methods can be found in Materials and methods) and asymptotic theory (Equation (9), black lines). Here, D=10-10 m2/s and v=2 µm/s with Cth/a⁢ρ=2/π⁢v. The initial signaling colony is of size ri=4⁢D/v (dashed vertical line). (C) Numerical wave speed as measured at t=100⁢D/v2 (markers) agrees well with theory (Equation (8), black line) as we independently vary D (circles) and Cth/a⁢ρ (diamonds) relative to the panel B values (blue circle and diamond). As predicted, v is indeed D-independent in this system.
Figure 3 Wave initiation dynamics. (A) Schematic demonstrating wave initiation. Cells within some initial signaling volume of radius ri begin signaling at some rate a. The signaling wave is initiated when the concentration at nearby cells exceeds the threshold concentration, Cth. (B) Cells near the initial signaling volume participate in the emission according an activation function, f⁢(c). For instance, in the case of Hill function activation, f⁢(c)=[1+(Cth/c)n]-1. C: Initiation times for Heaviside activation, in which f⁢(c)=Θ⁢[c-Cth]. Numerics (thick colored lines) and approximate asymptotic theory (Equations (10) and (11), thin black lines) of the initiation time’s dependence on ri for cells and diffusion in 1D or 2D (left) or cells in 2D with diffusion in 3D (right). For cells and diffusion in 1D, Equation (10) provides a good approximation in the limits v⁢ri/D≪1 and v⁢ri/D≫1. Similarly, for cells and diffusion in 2D, Equation (11) governs the large and small v⁢ri/D limits. In both of these cases, the wave always initiates, but the initiation time can be orders of magnitude larger than D/v2 if ri≪D/v. For cells and diffusion in 3D (right), signaling waves do not initiate for v⁢ri/D<3. Here again, the asymptotic theory Equation (12) is in good agreement with numerics.
Figure 4 Application to neutrophil swarming. (A) Schematic of the simple diffusion model. Here, cells on the target (within ri) signal distant neighbors by continuously emitting a single signaling molecule. If the neighboring cells have a chemotactic response, they migrate toward the target with some noise – that is, some non-zero angle θ with respect to the target. Otherwise, they move around with no sustained directionality. (B) Experimental data (color plot) reproduced from Reátegui et al., 2017 showing the information wave front in neutrophil swarming experiments. By tracking the neutrophils in space and time, they observe highly directed motion of the neutrophils towards the target (pink) starting around t=200 s. There is a clear boundary in space and time – the information wave front – between the regions where cells migrate toward the target (pink) and jostle around with no particular direction (white and light blue). While a relay theory (black line) is consistent with the convex shape of the information wave front, simple diffusive signaling by only the cells on the target (gray line) is not. The diffusion constants for both models is D=1.25×10-10 m2/s. The threshold concentrations for the relay and simple diffusion models are Cth/a⁢ρ≈3.66×105 s/m and 2.91×104 s/m, respectively. The parameters for the relay model are chosen to fit the wave front by eye while the simple diffusion model parameters are chosen to give the same signaling distance at t=500 s. (C) Gradients created by signaling relays (black) and simple diffusion (gray) models in panel B. The dashed vertical lines indicate the location of the information wave front. As time increases from left to right, the relay signaling motif gives an information wave that signals cells faster than simple diffusion in the long time limit. Cells within the wave front (to the left of the dashed lines that indicate the wave fronts) experience significantly larger gradients when the cells utilize a relay, which may facilitate efficient chemotaxis.
Appendix 3—figure 1 Relay concentration profiles. (A) Signaling molecule concentration profiles for cells in 2D and diffusion in 2D. Here, we assume continuous emission and no signaling molecule decay (left panel), continuous emission and signaling molecule decay (middle panel), or pulsed emission with no signaling molecule decay (right panel). The decay constant for the middle panel is γ=v2/4⁢D while the pulse width is τ=2⁢D/v in the right panel. For the left, center, and right panels, the threshold concentration is calculated according to Equations (24), (41), and (38), respectively. As discussed in the main text, the concentration profiles flatten (with respect to the profile generated by continuous emission without decay) inside the wave front once decay or pulsed emission is accounted for. (B) Signaling molecule concentration profiles for the same cases as in A, but with diffusion in 3D. Compared to the case of continuous emission without decay, the concentration profiles flatten (when accounting for decay) or have a local maximum (in the case of pulsed emission).
Appendix 3—figure 2 Wave speed with pulsed emission or decay. (A) Wave speed v for square pulse emission by cells in 2D with diffusion in 2D as a function of pulse width τ. We normalize v by the τ→∞ wave speed of (a⁢ρ⁢D/Cth)1/2. At τ=τc=h⁢Cth/a⁢ρ, the wave speed goes to zero, and for shorter pulses there is no wave-like solution. (B) Wave speed for continuous emission by cells in 2D with diffusion in 2D, accounting for signaling molecule decay at rate γ. For γ=tmin, 2D=1/2⁢τc, with τc as in panel A, the wave speed goes to zero and there are no wave-like solutions for larger decay rates. (C) Wave speed v as a function of pulse width τ for square pulse emission of a signaling molecule by cells in 2D with a 3D diffusive environment. At τ≈1.53⁢τc=1.53⁢D⁢(2⁢Cth/π⁢a⁢ρ)2 (vertical dashed line), there is a minimum wave speed of v≈0.6×(2⁢a⁢ρ/π⁢Cth) (horizontal dashed line), unlike with 2D diffusion. Importantly, τ≈1.53⁢τc is longer than the minimum initiation time of 4⁢τc/π. Thus, for values of τ between these two, cells can reach the threshold concentration but cannot propagate a traveling information wave. (D) Wave speed v as a function of decay rate γ for cells in 2D with diffusion in 3D. At γ=(π/4)2/τc, with τc as in C, the wave speed goes to zero.
Appendix 4—figure 1 Wave speed in a finite extracellular medium. (A) Method of images configuration for solving for concentrations in a finite-sized extracellular medium of thickness h. A point-like source (a cell) at z=0 is placed in an extracellular medium of thickness h. In order to calculate the concentration with the appropriate boundary conditions of ∂⁡c/∂⁡z at z=0,h, one need only add the contributions from ‘image cells’ at z=±2⁢j⁢h for j=1,2,…. (B) Universal curve (black line) showing non-dimensionalized wave speed (v⁢h/2⁢D) versus non-dimensionalized threshold concentration (π⁢Cth⁢D/a⁢ρ⁢h). In the limit of v⁢h/2⁢D≫1, we recover the familiar 3D scaling law of Equation (29) (blue line). In the limit of v⁢h/2⁢D≪1, we get the 2D scaling law of Equation (24) (red line).
Appendix 5—figure 1 Wave speeds and profiles with Hill function activation. (A) Numerical simulation of cells in 1D with diffusion in 1D and Hill function activation. The details of the simulation are described in Appendix 5: Asymptotic wave dynamics with Hill function activation. For n≥2, we see good agreement between the Heaviside theory (black lines) and the Hill function numerics (red dots). Snapshots are shown at t/τc=5,10,15,20,25 where τc=h2⁢Cth/a⁢ρ equals D/vΘ2, the characteristic time scale for Heaviside activation. Note that the x-axis is scaled by the characteristic length lc=h⁢(D⁢Cth/a⁢ρ)1/2, which is the length scale for Heaviside activation with a delta function source, D/vΘ. In the insets, we display the wave speed for the order-n Hill function, vn, compared to the Heaviside activation wave speed, vΘ=(a⁢ρ⁢D/h2⁢Cth)1/2. Our fit to the n=1 data gives vn=1≈1.93⁢vΘ, but we have shown in Appendix 5: Asymptotic wave dynamics with Hill function activation that vn=1=2⁢vΘ, meaning that the wave speeds in the insets are slight underestimates. (B) Numerical simulation of cells in 1D and diffusion in 2D with Hill function activation. For n≥2, the wave speed and concentration profiles (blue dots) agree well with the Heaviside theory (black lines). The theory plotted here assumes a delta-function-like source with respect to the extra diffusive dimension, as in Equation (25). The numerics, however, use a very narrow (H=lc/10) Gaussian source. Here, the characteristic length lc=π⁢h⁢D⁢Cth/2⁢a⁢ρ equals D/vΘ, the length scale for Heaviside activation. The x-axis is scaled by the same quantity. Snapshots are shown at t/τc=5,10,15,20,25 with the characteristic time τc=D⁢(π⁢h⁢Cth/2⁢a⁢ρ)2 equal to D/vΘ2, the time scale for Heaviside activation. In the insets, we display the wave speed for the order-n Hill function, vn, compared to the Heaviside activation wave speed, vΘ=2⁢a⁢ρ/π⁢Cth.
Appendix 6—figure 1 Wave dynamics with discrete sources. (A) Cartoon showing the discrete lattice model of cell signaling relays for cells in 1D. Here, a concentration wave propagates rightward at speed v through a group of cells spaced by a constant distance d. The concentration at x~=t=0 defines the threshold concentration in this discrete system and is a function of v,D, and d: Cth=Cth⁢(v,D,d). (B) Wave speed v compared to the continuous theory value of (a⁢D/d⁢h2⁢Cth)1/2 as a function of v⁢d/4⁢D for cells in 1D with diffusion in 1D. We can see that the wave speed approaches the continuous theory value for v⁢d/4⁢D≪1. (C) Same as B but for cells in 1D and diffusion in 2D.
Appendix 8—figure 1 Diffusive wave initiation with Hill function activation. (A) Wave initiation for cells in 1D with diffusion in 1D. Here, the initiation time tinit is normalized by the characteristic time scale τc=h2⁢Cth/a⁢ρ and the initial signaling colony size ri is normalized by the characteristic length scale lc=(h2⁢Cth⁢D/a⁢ρ)1/2. Yellow data points are individual simulations with dashed black guides to the eye connecting the yellow points. The solid yellow line corresponds to Heaviside-like activation and is reproduced from the main text. (B) Same as A, but for cells in 2D and diffusion in 2D. Here, τc=h⁢Cth/a⁢ρ and lc=(h⁢Cth⁢D/a⁢ρ)1/2. (C) Same as A, but for cells in 3D and diffusion in 3D. Here, τc=Cth/a⁢ρ and lc=(Cth⁢D/a⁢ρ)1/2. We note the value of ri=3⁢lc, the value below which waves fail to initiate for Heaviside-like activation.
Appendix 9—figure 1 Wavefront fitting of the simple diffusion and relay models. (A) Plot showing the information waves for various choices of diffusion constant (black) overlaid on the experimental chemotactic index data from Reátegui et al., 2017 (color plot). Here, we take ri=100 µm, although the target in the experiment is a smaller, oblong object. The size of the target has no effect on the convex shape of the information wave front. The black line is reproduced from the main text and has D=1.25×10-10 m2/s and asymptotic wave speed v≈1.7 µm/s (the threshold concentration is Cth/a⁢ρ=2/π⁢v). This diffusion constant is consistent with a small molecule like LTB4, and the resulting information wave dynamics can be made to fit the information wave front – both the initiation time of ≈200 s and the transient dynamics. Other choices of parameters (green: D=1.8×10-10 m2/s, v≈2.3 µm/s, navy: D=0.8×10-10 m2/s, v≈1.3 µm/s) give information wave fronts that are also roughly consistent with the experimental data. (B) However, with a much larger (D=10-9 m2/s, red) or smaller (D=10-11 m2/s, dashed) diffusion constant, an information wave with the correct initiation time does not have the correct transient dynamics. The wave speeds for these larger and smaller diffusion constants are, respectively, v≈11 µm/s and v≈0.3 µm/s. (C) Simple diffusion models of various diffusion constants overlaid atop data from Reategui et al. in which the the same swarming assay as in A/B (but with a slightly smaller 80 µm target) was utilized but with neutrophils whose LTB4 receptors (BLT1/2) had been blocked. Here, we see a 250 s offset before the propagation of a slow-moving diffusive front. A simple diffusive model does not capture this offset well, but does accurately capture the concave shape of the front (which contrasts with the convex shape of the front propagated by a relay). (D) Same as in C but with a 250 s offset in the theoretical curves. These curves fit the observed chemotactic index dynamics fairly well.
Appendix 10—figure 1 Comparison of relay and simple diffusion for cells and diffusion in 2D. All plots are for D=10-10 m2/s, v=2 µm/s, and ri=4⁢D/v. (A) Information fronts for relay (black) and simple diffusion (gray) models. The information wave travels like t for simple diffusion and v⁢t for the relay. The threshold concentration for the simple diffusion model is chosen such that the its information front intersects the relay model at t=20⁢D/v2. Thus, hCth,relay/aρ=25 s while h⁢Cth,diff./a⁢ρ≈0.25 s. (B) Snapshots of radial gradients at various times for both the relay (black) and simple diffusion (gray) models. The dashed vertical lines indicate the location of the wave fronts. For short times (left) the relay’s information wave front lags behind the simple diffusion model’s information front, though it later catches up (middle) and passes it (right). At all times, for cells just inside the wave front, the relay model creates gradients that are orders of magnitude larger than does simple diffusion.
Appendix 11—figure 1 Effects of chemotaxis on wavefront location and concentration profiles. (A) Information wave fronts for cell signaling relays with (navy) and without (black) chemotaxis. The information wave fronts are overlaid on the experimental chemotactic index data from (color plot) (Reátegui et al., 2017). The black curve is reproduced from the main text. Both models can account for the observed information wave fronts by fitting two parameters: the signaling molecule diffusivity, D and the threshold concentration, Cth. (B/C) Concentration profiles (B) and gradients (C) generated by the signaling relay models in A. The wave front is indicated in all panels by the dashed line, and the concentration profiles and gradients are plotted at the times such that the threshold concentration is equal to the concentration at the wave front. When one accounts for chemotaxis, the concentration profiles near the target steepen relative to models without chemotaxis.

Table 1

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