Flow-driven instabilities during pattern formation of Dictyostelium discoideum

In reaction–diffusion systems, an advective flow can induce unique emergent phenomena. One well known example is the differential flow induced chemical instability (DIFICI) that destabilizes an otherwise spatially homogeneous state of a system. DIFICI was originally predicted by Rovinsky and Menzinger [1, 2] and later confirmed experimentally for the Belousov–Zhabotinsky reaction [3, 4]. The basic idea behind DIFICI is that the reacting species flow at different speeds. This differential transport can initiate instabilities of the homogeneous steady state that induce propagating wave trains moving in the flow direction. Later, conditions for the emergence of three-dimensional spatial and spatiotemporal patterns due to DIFICI were derived by Malchow [5, 6] and demonstrated by patterns in a minimal phytoplankton–zooplankton model [7]. Instabilities in the homogeneous distribution can arise if phytoplankton and zooplankton move with different velocities, regardless of which one is faster. This mechanism of generating spatial structures is free from the restrictions of the Turing mechanism [8], which requires a large difference in diffusion coefficients of the two species involved. Accordingly one can expect DIFICI to be found widely in population dynamics [9] and biological morphogenesis [10].

Recently we conducted experiments to study the effect of a differential flow in a biological system, namely quasi one-dimensional colonies of the signaling amoebae Dictyostelium discoideum (D. d.) [11]. The aggregation of D. d. amoebae after nutrient deprivation is one of the best model systems for the study of spatial–temporal pattern formation at the multicellular level. Starved D.d. cells aggregate to form a multicellular structure by directional movement of cells in response to concentration changes of the chemoattractant cyclic adenosine monophosphate (cAMP) [12–14]. D. d. cells can be excitable [15, 16], and when stimulated with extracellular cAMP, they produce more cAMP, which in turn stimulates neighboring cells [17, 18]. This excitable response explains the wavelike pattern of cAMP, which propagates in the form of either concentric or spiral waves from the self-organized aggregation centers. In the natural environment, this aggregation occurs in forest soil and can be subject to water flow. In general, this medium into which the cAMP is released is a moving fluid while the cells are crawling on a solid substrate and are not transported by the flow.

In our recent experiments with chemotactically active populations of D.d. cells [11], a narrow straight microfluidic channel was used and cells were allowed to settle on the glass substrate before a laminar flow of buffer solution through the channel was started. This flow was not strong enough to detach the cells from the substrate and extracellular cAMP molecules were advected downstream. This differential transport in the extracellular medium induced macroscopic wave trains that developed spontaneously, propagated in the flow direction and had a unique period.

In this work, we have investigated the mechanisms of flow-driven waves using the well-established two and three component Martiel–Goldbeter (MG) model [19]. In the linear regime, our analytical calculations show that the convective transport of extracellular cAMP in a uniform field of signal-relaying cells leads to a flow-induced instability of the traveling-wave type. In the nonlinear regime, numerical simulations show propagating waves initiated by a convective instability, in agreement with the predictions of the linear analysis.

The kinetics of the cAMP production in well-stirred suspensions of D. d. cells are well described by the MG model [19], which is based on the relevant kinetic rate laws and the interaction between cAMP and its membrane receptor (see figure 1). In this model, extracellular cAMP binds to the receptor, and the receptor-cAMP complex activates adenylate cyclase which is an enzyme that produces intracellular cAMP from ATP. The newly synthesized cAMP is excreted into the extracellular medium where it binds to active receptors on the membrane and provides a feedback loop. A constant build up and accumulation of cAMP in the extracellular medium is limited by desensitization of the cAMP receptor upon prolonged exposure to high concentrations of cAMP. The receptor converts to an inactive form and loses its capacity to stimulate adenylate cyclase and synthesize more cAMP. Meanwhile extracellular phosphodiesterase which is supplied by the cells at a steady rate [20] degrades the cAMP in the system and in the absence of cAMP, the receptors slowly recover back to the active form. Note that the MG model is a mean-field type approximation and stochastic fluctuations are essentially neglected. The dynamics of the signaling system is given by the following set of coupled nonlinear ordinary differential equations

$$\begin{eqnarray}&&{{\partial }_{t}}\rho =-{{f}_{1}}(\gamma )\rho +{{f}_{2}}(\gamma )(1-\rho ),\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&{{\epsilon }_{2}}{{\partial }_{t}}\beta ={{s}_{1}}\Phi (\rho ,\gamma )-\beta ,\end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&{{\epsilon }_{1}}{{\partial }_{t}}\gamma ={{s}_{2}}\beta -\gamma ,\end{eqnarray} \tag{ 1c }$$

with

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {{f}_{1}}(\gamma ) & = & \frac{1+\kappa \gamma }{1+\gamma },\;\;\;\;\;\;\;\;{{f}_{2}}(\gamma )=\frac{{{L}_{1}}+\kappa {{L}_{2}}c\gamma }{1+c\gamma }, \\ \Phi (\rho ,\gamma ) & = & \frac{{{\lambda }_{1}}+{{Y}^{2}}}{{{\lambda }_{2}}+{{Y}^{2}}},\;\;\;\;\;\;\;\;Y=\frac{\rho \gamma }{1+\gamma }. \\ \end{array}\end{eqnarray*}$$

where t is equal to the real time multiplicated by k₁, ρ is the fraction of receptors in the active state, $\beta$ = ${{[{\rm cAMP}]}_{{\rm intracellular}}}/{{K}_{R}}$ and $\gamma \;=\;$ ${{[{\rm cAMP}]}_{{\rm extracellular}}}/{{K}_{R}}$ are dimensionless concentrations of intracellular and extracellular cAMP, respectively. Here ${{K}_{R}}={{10}^{-7}}\;{\rm M}$ is the dissociation constant of cAMP-receptor complex in an active state and ${{k}_{1}}=0.09\;{{{\rm min} }^{-1}}$ is the desensitization rate of active receptors. The dimensionless parameters are ${{\epsilon }_{1}}=\frac{{{k}_{1}}}{{{k}_{e}}}$, ${{\epsilon }_{2}}=\frac{{{k}_{1}}}{{{k}_{i}}+{{k}_{t}}}$, $\kappa =\frac{{{k}_{2}}}{{{k}_{1}}}$, ${{s}_{1}}=\frac{q\sigma }{{{k}_{t}}+{{k}_{i}}}\frac{\alpha }{1+\alpha }$, ${{s}_{2}}=\frac{{{k}_{t}}}{{{k}_{e}}h}$, ${{\lambda }_{1}}=\frac{\lambda \theta }{\epsilon }$, ${{\lambda }_{2}}=\frac{1+\alpha \theta }{\epsilon (1+\alpha )}$, ${{L}_{1}}=\frac{{{k}_{-1}}}{{{k}_{1}}}$, ${{L}_{2}}=\frac{{{k}_{-2}}}{{{k}_{2}}}$ and $c=\frac{{{K}_{R}}}{{{K}_{D}}}$. The values of model parameters are ${{k}_{2}}=1.665\;{{{\rm min} }^{-1}}$, ${{k}_{t}}=0.9\;{{{\rm min} }^{-1}}$, ${{k}_{e}}=2.5-12.5\;{{{\rm min} }^{-1}}$, ${{k}_{i}}=1.7\;{{{\rm min} }^{-1}}$, $\sigma =0.05-1.35\;{{{\rm min} }^{-1}}$, ${{L}_{1}}=10,{{L}_{2}}=0.005$, $\alpha =3$, $\theta =0.01$, $\lambda =0.01$, $\epsilon =1$, $h=5-100$, $c=10$, ${\rm and}$ $q=4000$. These parameter values are the experimental data listed in the last column of table II in [17], except for the values of k₁ = 0.09 min⁻¹ and k₂ = 1.665 min⁻¹ which are adjusted to obtain waves with a period of approximately 5 min. Here σ corresponds to the maximum activity of adenylate cyclase, which is the enzyme producing intracellular cAMP after being activated by an active cAMP-receptor complex, and k_e is the degradation rate of extracellular cAMP by the enzyme phosphodiestrase. Lastly, h is the ratio of the extracellular to intracellular volume which varies according to the cell density. For the definition of the other parameters see figure 1 and [19].

Zoom In Zoom Out Reset image size

Tyson and coworkers extended the MG model to a field of stationary signaling cells by adding spatial diffusion of cAMP through the extracellular medium [21]. The changes in both intracellular cAMP and the membrane receptors occur only by reactions described by equations (1a), (1b). However, in the presence of an external flow, the extracellular cAMP spreads out in the extracellular medium by both molecular diffusion and advection. The reaction–diffusion–advection form of the equation (1c) reads

$$\begin{eqnarray}&&{{\epsilon }_{1}}{{\partial }_{t}}\gamma ={{s}_{2}}\beta -\gamma +\epsilon _{1}^{2}{{\nabla }^{2}}\gamma +{{\epsilon }_{1}}\;v\nabla \gamma .\end{eqnarray} \tag{ 2 }$$

∇² and ∇ are the one-dimensional Laplacian and gradient operator for dimensionless spatial variable x defined by $x=\frac{{{k}_{1}}}{\sqrt{{{k}_{e}}D}}\times x^{\prime}$ from the dimensional space variable x'. $D=0.024\;{\rm m}{{{\rm m}}^{2}}{{{\rm min} }^{-1}}$ [21] is the diffusion coefficient of cAMP and v is the dimensionless flow velocity $v={{V}_{f}}/\sqrt{{{k}_{e}}D}$, where V_f is the dimensional velocity. Here we assume that the flow velocity V_f is well below the critical velocities needed to detach the cells from the substrate. Therefore the cells are not transported with the flow and only extracellular cAMP is advected downstream. Extracellular free phosphodiestrase is also subjected to the flow and therefore one would expect it to be spatially inhomogeneous or have insufficient concentration. Nevertheless, it is known that phosphodiestrase also occurs in membrane-bound form [22–24] and thus will be active also in the presence of an external flow.

To obtain a detailed understanding of flow-induced instabilities, we analyze the MG model both analytically and by numerical simulations. The analysis of the MG model is simplified if we eliminate the equation for β (intracellular cAMP) in favor of ρ (receptor) and γ (extracellular cAMP) [19, 21]. In the limit of ${{\epsilon }_{2}}\ll {{\epsilon }_{1}}$, this simplification is possible which gives $\beta \sim {{s}_{1}}\Phi (\rho ,\gamma )$. Then equations (1) in the presence of diffusion and external flow become

$$\begin{eqnarray}&&{{\partial }_{t}}\gamma ={{\epsilon }_{1}}{{\nabla }^{2}}\gamma +v\nabla \gamma +\frac{1}{{{\epsilon }_{1}}}[s\Phi (\rho ,\gamma )-\gamma ],\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&{{\partial }_{t}}\rho =-{{f}_{1}}(\gamma )\rho +{{f}_{2}}(\gamma )(1-\rho ),\end{eqnarray} \tag{ 3b }$$

where $s={{s}_{1}}{{s}_{2}}$.

We linearize equation (3) about the steady state $({{\gamma }_{0}},{{\rho }_{0}})$ writing

$$\begin{eqnarray}&&\gamma ={{\gamma }_{0}}+{{\gamma }_{1}},\;\;\;\;\;\;\;\;\rho ={{\rho }_{0}}+{{\rho }_{1}},\end{eqnarray} \tag{ 4 }$$

and assuming that ${{\gamma }_{1}}$ and ${{\rho }_{1}}$ are small. This gives

$$\begin{eqnarray}&&{{\partial }_{t}}{{\gamma }_{1}}={{a}_{11}}{{\gamma }_{1}}+{{a}_{12}}{{\rho }_{1}}+{{\epsilon }_{1}}{{\nabla }^{2}}{{\gamma }_{1}}+v\frac{\partial }{{{\partial }_{x}}}{{\gamma }_{1}}\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\rho }_{1}}={{a}_{21}}{{\gamma }_{1}}+{{a}_{22}}{{\rho }_{1}},\end{eqnarray} \tag{ 5b }$$

where

$$\begin{eqnarray}&&{{a}_{11}}=\frac{s{{N}_{0}}-1}{{{\epsilon }_{1}}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{a}_{12}}=\frac{s{{M}_{0}}}{{{\epsilon }_{1}}},\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{{a}_{22}}=-{{f}_{1}}({{\gamma }_{0}})-{{f}_{2}}({{\gamma }_{0}}),\;\;\;\;\;\;\;\;\;\;\;\;{{a}_{21}}=(1-{{\rho }_{0}}){{A}_{0}}-{{B}_{0}}{{\rho }_{0}},\end{eqnarray} \tag{ 6b }$$

and

$$\begin{eqnarray}&&{{M}_{0}}=\frac{2{{\rho }_{0}}\gamma _{0}^{2}({{\lambda }_{2}}-{{\lambda }_{1}})}{{{(1+{{\gamma }_{0}})}^{2}}{{({{\lambda }_{2}}+Y_{0}^{2})}^{2}}},\;\;\;\;\;{{N}_{0}}=\frac{2\rho _{0}^{2}{{\gamma }_{0}}({{\lambda }_{2}}-{{\lambda }_{1}})}{{{(1+{{\gamma }_{0}})}^{3}}{{({{\lambda }_{2}}+Y_{0}^{2})}^{2}}},\;\;\;\;\;{{Y}_{0}}=\frac{{{\rho }_{0}}{{\gamma }_{0}}}{1+{{\gamma }_{0}}},\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{{A}_{0}}=\frac{(\kappa {{L}_{2}}-{{L}_{1}})c}{{{(1+c{{\gamma }_{0}})}^{2}}},\;\;\;\;{{B}_{0}}=\frac{\kappa -1}{{{(1+{{\gamma }_{0}})}^{2}}}.\end{eqnarray} \tag{ 7b }$$

Note that the stability of the system in the absence of diffusion and flow requires that the trace of system (5) denoted by $T={{a}_{11}}+{{a}_{22}}$ and its determinant $\Delta ={{a}_{11}}{{a}_{22}}-{{a}_{12}}{{a}_{21}}$ should satisfy the following conditions

$$\begin{eqnarray}&&T\lt 0,\;\;\;\;\;\;\;\Delta \gt 0.\end{eqnarray} \tag{ 8 }$$

We represent the perturbation ${{\gamma }_{1}}$ and ${{\rho }_{1}}$ by the spatial Fourier expansion

$$\begin{eqnarray}&&({{\gamma }_{1}},{{\rho }_{1}})=\frac{1}{2\pi }\int _{-\infty }^{\infty }({{\gamma }_{k}},{{\rho }_{k}}){{{\rm e}}^{{\rm i}kx+\omega (k)t}}{\rm d}k\end{eqnarray} \tag{ 9 }$$

where ω satisfies the dispersion relation

$$\begin{eqnarray}&&{{\omega }^{2}}+\omega (-T+{{\epsilon }_{1}}{{k}^{2}}+{\rm i}vk)+\Delta -{{a}_{22}}({{\epsilon }_{1}}{{k}^{2}}+{\rm i}vk)=0.\end{eqnarray} \tag{ 10 }$$

Note that diffusion alone cannot destabilize the system, because for v = 0, conditions (8) and ${{a}_{22}}=-{{f}_{1}}({{\gamma }_{0}})-{{f}_{2}}({{\gamma }_{0}})\lt 0$ imply that $\operatorname{Re}(\omega )$ is negative for all wavenumbers k. However, a reaction–diffusion-flow system may be unstable for nonzero v. To investigate this flow-driven instability, we first plot $\operatorname{Re}(\omega )$ as a function of wavenumber for different flow velocities ${{V}_{f}}({\rm mm}\;{{{\rm min} }^{-1}})=v\sqrt{{{k}_{e}}D}$. If ${{a}_{11}}\gt 0$, then $\operatorname{Re}(\omega )$ becomes positive at sufficiently large k and the system becomes unstable to short-wavelength perturbations. However, the diffusion term produces a short-wavelength cutoff in the dispersion curves as shown in figure 2(a). This entails the appearance of a threshold flow velocity v_min, below which the homogeneous steady state is always stable. To obtain the critical velocity, we look at the neutral curve in the $v-{{k}^{2}}$ plane on which $\operatorname{Re}(\omega )=0$. This curve is obtained from the dispersion relation (10) and is given by

$$\begin{eqnarray}&&{{v}^{2}}=\frac{(-\frac{\Delta }{{{a}_{22}}}+{{\epsilon }_{1}}{{k}^{2}}){{({{\epsilon }_{1}}{{k}^{2}}-T)}^{2}}}{{{k}^{2}}({{a}_{11}}-{{\epsilon }_{1}}{{k}^{2}})}.\end{eqnarray} \tag{ 11 }$$

Note that equation (11) needs the additional condition ${{a}_{11}}\gt 0$ to be defined on the interval $0\lt {{\epsilon }_{1}}{{k}^{2}}\lt {{a}_{11}}$. A typical shape of the neutral curve is presented in figure 2(b). It is possible to obtain a critical velocity as well as a critical wave number from the conditions $\operatorname{Re}[\omega ({{k}_{c}})]=0$ and ${\rm d}[\operatorname{Re}\omega ({{k}_{c}})]/{\rm d}{{k}^{2}}\ \;=\ \;0$. A rough estimate gives ${{k}_{c}}\sim {{(|{{a}_{11}}+{{a}_{22}}|/{{\epsilon }_{1}})}^{1/2}}$ and ${{v}_{{\rm min} }}\sim {{[({{a}_{11}}+{{a}_{22}})\Delta {{\epsilon }_{1}}/{{a}_{11}}{{a}_{22}}]}^{1/2}}$. Note that ${{v}_{{\rm min} }}$ depends on all the parameters in the model. Exemplarily, figure 2(c) shows that the threshold velocity ${{V}_{{\rm min} }}$ is smaller for low cell densities (high h values). The comparison between linear stability analysis and numerical simulations is presented in figure 2(d), where the values of k with maximum growth rate are plotted for different flow velocities. Interestingly, there is a good agreement between simulations and linear stability analysis especially at large and small values of flow rates.

Zoom In Zoom Out Reset image size

As shown by the bifurcation diagram in figure 3, the parameters σ, k_e and h are among the important factors governing the dynamics of the MG model (equations (3)). The region of possible convective instability is labeled CU in the figure where the system can sustain convective instabilities for flow velocities above the critical ${{v}_{{\rm min} }}$. The results in figure 3(a) predict that convective instabilities can occur for both large and small values of h which correspond to lower and higher cell densities, respectively. However, the 'window' for convective instabilities is slightly narrower for larger values of h (low cell density). More importantly, figure 3(b) predicts that, for fixed values of the other parameters, there is a finite range of degradation rate k_e just before the absolute unstable (AU) regime, where convective instabilities can occur. This window increases slightly as σ increases (i.e. higher activity of the enzyme adenylate cyclase).

Zoom In Zoom Out Reset image size

We solve equations (3) numerically using the Merson modification of the Runge–Kutta method with automatic regulation of the time step. The spatial domain consists of 6000 grid points in the x direction with a mesh size of $\Delta x=0.025$ mm. The flow is from the left to the right and a no-flux boundary condition is used on the right boundary that corresponds to the outlet of the channel. To have a realistic description of the experiment at the left side of the channel [11], we assume that over the first 2.5 mm of the channel length (corresponding to the inlet of the microfluidic channel), there are no cells and the active production of cAMP drops to zero. cAMP diffuses back from the channel to the inlet area which can then be advected downstream. This boundary condition at the inlet has a significant effect on the observed pattern as will be discussed later.

To show the development of a convective instability along the channel, we select k_e = 3 min⁻¹ and $\sigma =0.45$ min⁻¹ to be in the CU regime of the phase diagram, but close to transition to AU regime (point $*2$ in figure 3(b)). The flow-driven waves are sustained at flow velocities greater than the critical velocity v_min. In the inlet region, we assume that initially the concentration of cAMP is zero. At t = 0, the spatially uniform steady state of the system is perturbed by slightly increasing the value of γ over a small region in the close neighborhood of the inlet (figure 4(a)). No additional perturbations are applied during the course of the integration. The perturbation applied at t = 0 generates the first wave which propagates along the channel (figure 4(b)). The second wave appears in the downstream of the location where the initial perturbation was applied (figure 4(c)). Note that the new pulse emerge from the back of the already existing pulse, as shown in figures 4(c) and (f). The effect of the refractory tail of the 'mother' pulse is to reduce the amplitude of the 'daughter' wave. As the two wave fronts become more separated, the 'daughter' wave grows to form a fully developed new 'mother' wave. An important feature to note is that each successive wave is formed further downstream than the previous one and the old waves die out as they propagate through the channel (see figure 4). The growth of new wave fronts and decay of the old ones is clearly visible in the corresponding space-time plot presented in figure 5(b). The traveling waves are advected downstream and eventually out of the system with time. The waves initially propagate with velocity of 2.20 mm min⁻¹, which is larger than the imposed flow velocity and as the wave amplitude decays, they slow down to 1.75 mm min⁻¹. Interestingly, the period of the flow-induced waves decrease from 5 to 2.4 min. Similar behavior is observed in figure 5(a) where the parameters are selected to be deep inside the CU regime (point $*1$ in figure 3(b)). Here convective instability waves have smaller amplitudes and the wave package extends over a smaller region of the channel. The bottom two plots in figures 5(c),(d) correspond to simulations in the AU and excitable regime, respectively. The flow driven waves in figure 5(c) are also traveling downstream and interact with uniform bulk oscillations at the lower part of the channel and sweep them away. Note that the boundary condition at the upstream of the channel (zero cAMP concentration) suppresses the bulk oscillations at the upper part of the channel. Similarly, the wave velocity varies from 2.5 to 1.95 mm min⁻¹ before interaction with bulk oscillations and the wave period varies from 5 to 2.5 min. In the excitable regime, if the perturbation amplitude at t = 0 is larger than a critical value, then it is amplified and one single excitation wave propagates along the channel (figure 5(d)), otherwise the perturbation dies out.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows the results obtained at lower flow velocities, namely V_f = 1 mm min⁻¹. In figure 6(a), the flow velocity is below the critical velocity ${{v}_{c}}\sim 1.30$ mm min⁻¹ and therefore convective instability waves are not observed. However, the cell-free boundary region at the inlet drives a new kind of instability in the area close to the boundary. cAMP secreted by the cells diffuse back to the inlet and at the same time can be advected downstream by the flow. Since there are no cells in the inlet to produce cAMP, the points in the vicinity of the inlet lose cAMP due to advection or diffusion and gain only a little from the upstream of the channel (inlet). In other words, there is a large negative flux of cAMP in the neighborhood close to the inlet, which can be considered as a sink. Under some flow conditions, this sink drives the system from a stable to an oscillatory regime. To obtain a better understanding of sink driven oscillations (SDOs), we performed linear stability analysis of equation (3) for an isolated point with a constant negative flux of cAMP due to advection and diffusion. We also plot the nullclines of equation (3) with zero flux (i.e. no diffusion and no advection) and with negative constant flux of cAMP. Negative constant flux of cAMP shifts one of the nullclines upwards (equation (3a)), but does not change the second nullcline (equation (3b)). Linear stability analysis predicts that for negative flux smaller than a critical value f_c, the new fixed point is stable (figures 7(a)–(c)). However, it becomes unstable for a negative flux larger than f_c and the system is attracted to a limit cycle (see figure 7(d)). This negative flux can be considered as a driving force for SDOs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we run simulations for an extended system (equation (3)) with advection and diffusion for different flow velocities (figure 8). The results are shown for a point located very close to the inlet. At small flow velocities below 0.275 mm min⁻¹, the new fixed point is stable and oscillations in cAMP concentration are damped. It is important to emphasize that e.g. in figure 8(a), at the end the system has a steady negative flux of value −15.3. This negative flux is supposed to destabilize the system, in accordance with linear stability analysis of an isolated point. Interestingly, advection and diffusion stabilize possible oscillations and the system is attracted by the new fixed point. As we increase the flow velocity, at V = 0.275 mm min⁻¹, the shifted nullcline is tangent to the zero-flux nullcline and the system is at the onset of instability (figure 8(b)). Above this critical velocity there is no intersection point anymore, except at very small values of cAMP concentration and a stable limit cycle appears which leads to oscillations in cAMP concentration (see figure 8(c)). In other words, close to the inlet, there is a competition between the destabilizing effect of the inlet and stabilizing effect of diffusion and advection. However, there is an upper limit of flow velocity above which the limit cycle disappears and the system is attracted by the fixed point with small concentration of cAMP (figure 8(d)). For the set of parameters in figure 6(a), SDOs are observed for flow velocities larger than 0.2 mm min⁻¹ but smaller than 1.3 mm min⁻¹. The amplitude of these SDOs decay as the distance to the inlet region increases. Far away from the inlet, the net flux of cAMP is very close to zero and the new fixed point is stable.

Zoom In Zoom Out Reset image size

In figure 6(b), the SDOs are combined with flow driven waves, where the flow velocity is well above the critical velocity ${{v}_{c}}\sim 0.12$ mm min⁻¹. Similarly, a combination of SDOs and flow-driven waves is observed in figure 6(c), which is very similar to the experimental observations in [11]. In the experiments, the channel length was only 30 mm. Thus figure 6(c) predicts that for experiments with much longer channels, uniform bulk oscillations may be observed at the downstream of the channel.

In this section, we consider the full set of equations with diffusion and advection (equations (1a), (1b) and (2)). If we assume that the system has a stable steady state $({{\gamma }_{0}},{{\rho }_{0}},{{\beta }_{0}})$, then the evolution of a small perturbation is governed by the linearized equations

$$\begin{eqnarray}&&{{\partial }_{t}}{{\gamma }_{1}}={{a}_{11}}{{\gamma }_{1}}+{{a}_{12}}{{\rho }_{1}}+{{a}_{13}}{{\beta }_{1}}+{{\epsilon }_{1}}{{\nabla }^{2}}{{\gamma }_{1}}-v{{\partial }_{x}}{{\gamma }_{1}},\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\rho }_{1}}={{a}_{21}}{{\gamma }_{1}}+{{a}_{22}}{{\rho }_{1}}+{{a}_{23}}{{\beta }_{1}},\end{eqnarray} \tag{ 12b }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\beta }_{1}}={{a}_{31}}{{\gamma }_{1}}+{{a}_{32}}{{\rho }_{1}}+{{a}_{33}}{{\beta }_{1}},\end{eqnarray} \tag{ 12c }$$

where

$$\begin{eqnarray}&&\gamma ={{\gamma }_{0}}+{{\gamma }_{1}},\;\;\;\;\;\;\;\;\rho ={{\rho }_{0}}+{{\rho }_{1}},\;\;\;\;\;\;\;\;\beta ={{\beta }_{0}}+{{\beta }_{1}},\end{eqnarray} \tag{ 13 }$$

where ${{\gamma }_{1}},{{\rho }_{1}},{{\beta }_{1}}$ are small and a_ij are the elements of the Jacobian matrix

$$\begin{eqnarray}&&{{a}_{11}}=\frac{-1}{{{\epsilon }_{1}}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{a}_{12}}=0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{a}_{13}}=\frac{{{s}_{2}}}{{{\epsilon }_{1}}},\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&{{a}_{21}}=(1-{{\rho }_{0}}){{A}_{0}}-{{B}_{0}}{{\rho }_{0}},\;\;\;\;\;\;\;\;\;\;\;\;{{a}_{22}}=-{{f}_{1}}({{\gamma }_{0}})-{{f}_{2}}({{\gamma }_{0}}),\;\;\;\;\;\;\;\;\;\;\;\;{{a}_{23}}=0,\end{eqnarray} \tag{ 14b }$$

$$\begin{eqnarray}&&{{a}_{31}}=\frac{{{s}_{1}}N}{{{\epsilon }_{2}}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{a}_{32}}=\frac{{{s}_{1}}M}{{{\epsilon }_{2}}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{a}_{33}}=\frac{-1}{{{\epsilon }_{2}}},\end{eqnarray} \tag{ 14c }$$

with M₀, N₀, A₀ and B₀ defined in equation (7).

We first consider the system without diffusion and convection and introduce the following definitions: $T={{a}_{11}}+{{a}_{22}}+{{a}_{33}}$ is the trace of the system, Δ is its determinant, ${{T}^{ij}}={{a}_{ii}}+{{a}_{jj}}\;(i\lt j)$ is the trace of the subsystem formed by $\gamma ,\rho \;(\beta =0)$ or $\gamma ,\beta \;(\rho =0)$ or $\rho ,\beta \;(\gamma =0)$, and ${{\Delta }_{ij}}={{a}_{ii}}{{a}_{jj}}-{{a}_{ij}}{{a}_{ji}}\;(i\lt j)$ is the determinant of the same subsystem, while $\Sigma ={{\Delta }_{12}}+{{\Delta }_{13}}+{{\Delta }_{23}}$. Based on the Hurwitz criterion [25], the stability of the entire three variable system requires that

$$\begin{eqnarray}&&T\lt 0,\;\;\;\;\;\;\;\Delta \lt 0,\;\;\;\;\;\;{\rm and}\;\;\;\;\;\;\;\;\;\;T\Sigma -\Delta \lt 0.\end{eqnarray} \tag{ 15 }$$

Let us now consider the system (12) with flow but without the diffusion term. We represent the perturbations $({{\gamma }_{1}},{{\rho }_{1}},{{\beta }_{1}})$ by the spatial Fourier transform $({{\gamma }_{1}},{{\rho }_{1}},{{\beta }_{1}})=\frac{1}{2\pi }\int _{-\infty }^{\infty }({{\gamma }_{k}},{{\rho }_{k}},{{\beta }_{k}}){{{\rm e}}^{{\rm i}kx+\omega (k)t}}{\rm d}k$. The equations for the Fourier components then become

$$\begin{eqnarray}&&{{\partial }_{t}}{{\gamma }_{k}}={{a}_{11}}{{\gamma }_{k}}+{{a}_{12}}{{\rho }_{k}}+{{a}_{13}}{{\beta }_{k}}-v{\rm i}k{{\gamma }_{k}},\end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\rho }_{k}}={{a}_{21}}{{\gamma }_{k}}+{{a}_{22}}{{\rho }_{k}}+{{a}_{23}}{{\beta }_{k}},\end{eqnarray} \tag{ 16b }$$

$$\begin{eqnarray}&&{{\partial }_{t}}{{\beta }_{k}}={{a}_{31}}{{\gamma }_{k}}+{{a}_{32}}{{\rho }_{k}}+{{a}_{33}}{{\beta }_{k}}.\end{eqnarray} \tag{ 16c }$$

We form the characteristic polynomial $P({\rm i}\omega )={{P}_{R}}(\omega )+{\rm i}{{P}_{I}}(\omega )$, where ω, P_R, and P_I are all real. The real and imaginary parts of $P({\rm i}\omega )$ are

$$\begin{eqnarray}&&{{P}_{R}}(\omega )=T{{\omega }^{2}}+vk{{T}_{23}}\omega -\Delta ,\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{{P}_{I}}(\omega )=-{{\omega }^{3}}-vk{{\omega }^{2}}+\Sigma \omega +kv{{\Delta }_{23}}.\end{eqnarray} \tag{ 17b }$$

For convective instability to occur, the following conditions should be fulfilled [2]

$$\begin{eqnarray}&&b/{{a}_{11}}\gt 0\;\;\;\;\;{\rm and}\;\;\;\;{{b}^{2}}-4{{a}_{11}}\Delta {{\Delta }_{23}}\gt 0.\end{eqnarray} \tag{ 18 }$$

where $b=\Delta +{{a}_{11}}{{\Delta }_{23}}-{{T}_{23}}({{\Delta }_{12}}+{{\Delta }_{13}})$.

Including the diffusion of extracellular cAMP in addition to the convection term generates a short wavelength cutoff in the dispersion relation and leads to the appearance of a critical velocity, similar to the two-component case. We used the criteria in equation (18) to explore the regions of convective instabilities in the ${{k}_{e}}-h$ and $\sigma -{{k}_{e}}$ parametric planes (see figure 9). Note that again, similar to the two-component MG model, regions of CU occur right before the absolute instability (AU) but in a narrower region (figure 9(b)).

Zoom In Zoom Out Reset image size

We extended the numerical simulations described in section 3 to the three-component MG model and performed simulations for the set of parameters inside the CU regime, namely ${{k}_{e}}=2.5\;{{{\rm min} }^{-1}}$ and $\sigma =0.34\;{{{\rm min} }^{-1}}$. This point is marked by the ⋆ symbol in figure 9(b) and is located close to the transition line but in the convectively unstable region of the phase diagram. Note that, in the absence of flow (and diffusion), the homogeneous steady state is stable to small perturbations; however an external flow above the threshold velocity v_min may lead to convective instability. The corresponding wave pattern is illustrated in figures 10(a), (c). Remarkably, the flow driven instability waves of cAMP need more time to be fully developed. Since the system is very close to the transition line to the stable region, the wave fronts have a lower concentration of cAMP. However, deep inside the CU regime and close to the lower transition line in figure 9(b), the instability waves appear earlier in time and have higher concentrations of cAMP (see figures 10(b), (d)). This point is denoted by + in figure 9(b) and corresponds to parameter values ${{k}_{e}}=2.5\;{{{\rm min} }^{-1}}$ and $\sigma =0.24\;{{{\rm min} }^{-1}}$.

Zoom In Zoom Out Reset image size

Our linear stability analysis and numerical simulations of the MG model have shown that an advective flow of extracellular cAMP through a uniform field of signaling D.d. cells can induce a pattern forming instability. We found that, in the convectively unstable regime, small perturbations in the concentration of cAMP grow in the course of time and generate a wave train moving only in the flow direction. The flow-driven instability presents a non-oscillatory state without fluid flow that develops a pattern if the flow velocity increases beyond a critical threshold. We have repeated our numerical simulations for different flow rates and measured the wave propagation velocity as well as the wavelength. The calculated wave velocity and wavelength for different flow rates follow a linear relationship which is in agreement with experimental observations in [11]. However, the convective nature of this instability which is recognized in the numerical simulations as well as in linear stability analysis, is not directly observed in the experiments. Moreover, there are important experimental features that are different from numerical simulations. In the experiments, the propagating flow-driven waves appear to fill the whole length of microfluidic channel. This is different from simulations where each new wave front is formed further downstream than the previous one. Perhaps this is due to small continuous inlet-perturbations which exist throughout the experiment and can be considerably amplified by convective instabilities. These time-dependent inlet perturbations are not included in the numerical simulations where only a single perturbation was applied as an initial condition. Another possible explanation is the effect of the boundary which leads to a new kind of instability close to the inlet. As discussed in great detail in section 3, SDOs initiate the wave patterns in the vicinity of the inlet which is more similar to our experimental observations in [11]. This suggests that wave patterns observed in the experiment, could be a combination of flow-driven waves and SDOs. In addition, there is a significant difference between the experimental system which is a living system under development and the system studied numerically. In the numerical simulations, the point selected in the convectively unstable regime of the phase diagram is stable when the external flow is switched off and always remains at the same place in the phase diagram. However, D.d. amoebae represent a system, which is under continuous development and in the absence of external flow, they synchronize and show formation and propagation of spiral waves. In other words, in the course of development, the living system can have a developmental path in the phase diagram as has been proposed in [26].

There are relevant biological questions related to the influence of flow-driven waves on aggregation pattern of D.d. cells. These waves are observed in the initial stages of cell signaling, where the cell population is assumed to have a uniform distribution. However as the cells start to aggregate, variations in the cell density become significant. Future studies should attempt to include chemotactic cell movement in the simulations to study the distribution of cell population in the presence of flow-driven waves. Our results may also provide a starting point for the study of pattern formation and locomotion-driven instabilities in other biological systems.

We are grateful to Alan Pumir for discussions and comments. A G was supported by Dorothea-Schlözer scholarship from Georg-August University of Göttingen, Germany. O St acknowledges support provided by the A v Humboldt Foundation and the National Science Foundation (CHE-1213259).