Nonlinear self-adapting wave patterns

One of the ways in which a non-equilibrium system can lead to pattern formation is via a traveling wave bifurcation [1]. In such a system, the uniform state becomes unstable to modes at a finite wave vector k and a finite frequency ω leading to a variety of phenomena involving nonlinear traveling wave states. This scenario has proven relevant for processes ranging from binary fluid convection [2, 3] to electro-hydrodynamics in liquid crystals [4] to the sloshing of Min proteins in bacteria [5–14].

To date, this instability has been viewed as an oscillatory analog of the familiar Turing instability. Such an analogy implies that as a function of some control parameter, the instability sets in at a critical threshold value of the parameter at some given $k\ne 0$. Furthermore, above threshold, there is a band of unstable wavevectors stretching from $0\lt {k}_{\min }\lt k\lt {k}_{\max }$. Here we show that there exists a previously unconsidered possibility, namely that for some systems with a particular symmetry, k_min equals zero. This dramatically changes the nature of the nonlinear patterns that form, as there is no predetermined length scale for the emergent structure; instead, the waves are able to self-adapt to the size of the physical system. As we will discuss, this is the traveling wave analog of what happens in viscous fingering [15, 16] where the static bifurcation extends to ${k}_{\min }=0$. Models for the aforementioned Min dynamics offer a specific realization of this new paradigm. Moreover, the self-adaptation provides a mechanism whereby the dynamical pattern can maintain a one node form as the cell expands during growth.

We start with the model of [11, 17] for the Min system. There are two proteins, MinD and MinE, each of which can be on the membrane (m) or in the cytosol (c). The two proteins can reversibly desorb and adsorb, and adsorption of MinE involves it directly binding to an already membrane resident MinD. Additional nonlinearities emerge from the assumed cooperativity in the desorption rates. Adding in diffusion in the compartments leads in one spatial dimensions to the 4 coupled pde's

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{D}}{\partial t}={R}_{{de}\to D+E}-{R}_{D\to d}+{D}_{D}\displaystyle \frac{{\partial }^{2}{c}_{D}}{\partial {x}^{2}},\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{E}}{\partial t}={R}_{{de}\to D+E}-{R}_{d+E\to {de}}+{D}_{E}\displaystyle \frac{{\partial }^{2}{c}_{E}}{\partial {x}^{2}},\end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{d}}{\partial t}={R}_{D\to d}-{R}_{d+E\to {de}}+{D}_{d}\displaystyle \frac{{\partial }^{2}{c}_{d}}{\partial {x}^{2}},\end{eqnarray} \tag{ 1c }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{{de}}}{\partial t}={R}_{d+E\to {de}}-{R}_{{de}\to D+E}+{D}_{{de}}\displaystyle \frac{{\partial }^{2}{c}_{{de}}}{\partial {x}^{2}},\end{eqnarray} \tag{ 1d }$$

where c_D and c_E are the cytosol concentrations of MinD and MinE, c_d is the concentration of MinD on the membrane and c_de is the concentration of the MinD/MinE complex on the membrane and the rates are

$$\begin{eqnarray}\begin{array}{rcl}{R}_{D\to d} & = & ({\lambda }_{D}+{\lambda }_{{dD}}{c}_{d}){c}_{D},\\ {R}_{{de}\to D+E} & = & {\lambda }_{{de}}{c}_{{de}},\\ {R}_{d+E\to {de}} & = & ({\lambda }_{E}+{\lambda }_{{eE}}{c}_{{de}}^{2}){c}_{d}{c}_{E}\mathrm{.}\end{array}\end{eqnarray} \tag{ 2 }$$

Of critical importance to the physics of this system are two features. First, the diffusion constants of the membrane-bound species are orders of magnitude smaller than their values for the same protein in the cytosol. We will see that this property will guarantee the existence of an instability which drives the pattern formation and also will enable a simplification of the dynamics, at least for small system size. Second, the model contains two global conservation laws; the total number of both MinD and MinE proteins in the system are unchanged by each of the reactions, so that they are determined solely by the initial conditions. To see what this implies, we imagine we have found a uniform steady-state solution of the model and calculate the fate of a small perturbation around this solution, so that

$$\begin{eqnarray*}&&{c}_{D}={c}_{D}^{(0)}+{\delta }_{D}{{\rm{e}}}^{{\rm{i}}{k}{x}+{\rm{\Omega }}t}\end{eqnarray*}$$

with analogous expressions for the other concentrations. This leads in a standard manner to a 4 × 4 homogeneous linear system, with the growth rates Ω being given by the eigenvalues of the k-dependent stability matrix. The results of such a calculation for the most unstable mode, where for the moment we ignore the quantization of k induced by the boundaries, are presented in figure 1. We see that for $k\lt {k}_{* }$, there is a complex unstable mode, whose (complex) growth rate goes to 0 in the limit $k\to 0$. We shall now show that this is a direct result of the global conservation laws. At k = 0, due to these laws, the sums of equations 1(a, b, d) and equations 1(b, d) are identically zero, as shifts in the overall levels of either the conserved MinD or MinE proteins due to the perturbation are left unchanged by the dynamics. Hence ${\rm{\Omega }}=0$ is a doubly degenerate eigenvalue for k = 0, with left eigenvectors ${\hat{\phi }}^{(1)}=(1,0,1,1)$ and ${\hat{\phi }}^{(2)}=(0,1,0,1)$. We can then calculate the shift in ${\rm{\Omega }}(k)$ for these two modes for small k using what amounts to degenerate perturbation theory in quantum mechanics, generalized to non-Hermitian matrices. To leading order, we computing the projection of the perturbation, i.e., the diffusion terms in the stability matrix, namely ${ \mathcal D }={\rm{diag}}(-{D}_{D},-{D}_{E},-{D}_{d},-{D}_{{de}}){k}^{2}$ onto the basis of the degenerate $2x2$ subspace to obtain the reduced stability matrix

$$\begin{eqnarray*}&&{{ \mathcal S }}_{{ij}}={\hat{\phi }}^{(i)}{ \mathcal D }{\phi }^{(j)},\end{eqnarray*}$$

where the ${\phi }^{(j)}$ are the corresponding right eigenvectors of the two k = 0 zero modes, which satisfy the orthonomality condition ${\hat{\phi }}^{(i)}{\phi }^{(j)}={\delta }_{{ij}}$. We then have to diagonalize ${{ \mathcal S }}_{{ij}}$ to find ${\rm{\Omega }}(k)$. For the parameter set used in [11], the steady-state is ${c}_{D}^{0}=71.77$, ${c}_{E}^{0}=76.39$, ${c}_{d}^{0}=604.62$, ${c}_{{de}}^{0}=323.61$ and this gives rise to

$$\begin{eqnarray*}{\phi }^{(1)} & = & \left(\begin{array}{c}-0.0767\\ -0.1494\\ 0.9273\\ 0.1494\end{array}\right)\qquad {\phi }^{(2)}=\left(\begin{array}{c}0.4022\\ -0.0422\\ -1.4444\\ 1.0422\end{array}\right)\\ -{ \mathcal S }/{k}^{2} & = & {D}_{c}\left[\begin{array}{cc}-0.077 & 0.402\\ -0.149 & -0.042\end{array}\right]+{D}_{m}\left[\begin{array}{cc}1.077 & -0.402\\ 0.149 & 1.042\end{array}\right]\end{eqnarray*}$$

for pair-wise equal cytosol (D_c) and membrane (D_m) diffusivities. This matrix always has a pair of complex eigenvalues, as

$$\begin{eqnarray}&&{{\rm{tr}}}^{2}(S)-4{\rm{\det }}(S)=-0.239{k}^{4}{({D}_{{\rm{c}}}-{D}_{{\rm{m}}})}^{2}\lt 0.\end{eqnarray} \tag{ 3 }$$

This complex pair is unstable as long as

$$\begin{eqnarray}&&{\rm{tr}}(S)={k}^{2}(-0.119{D}_{{\rm{c}}}+2.119{D}_{{\rm{m}}})\lt 0\Rightarrow \displaystyle \frac{{D}_{{\rm{c}}}}{{D}_{{\rm{m}}}}\gt 17.8\end{eqnarray} \tag{ 4 }$$

which is easily satisfied by the biophysical parameters, as we saw in figure 1. In general, the initial rise of ${\rm{Re}}{\rm{\Omega }}$ with k depends on the relatively large cytosol diffusion constants whereas the value of k at which the system restabilizes depends on the small membranal ones. In the limit of very large diffusion constant ratio between the cytosol and membranal fields and for equal membranal diffusivities, D_m, the spectrum approaches (for k strictly non-zero) the simple form ${{\rm{\Omega }}}_{0}-{D}_{{\rm{m}}}{k}^{2}$ for complex ${{\rm{\Omega }}}_{0}$ with a positive real part.

Zoom In Zoom Out Reset image size

This stability structure presents a new twist on what happens in pattern forming systems such as viscous fingering and dendritic crystal growth [15, 16]. There, translation invariance of the base system guarantees a single zero k = 0 eigenvalue which gives rise to a real-mode instability for $0\lt k\lt {k}^{* }$. A related idea has arisen in the context of cellular processes that have one chemical component being exchanged between different compartments but is globally conserved [18]. In our system the existence of two zero modes and of course the non-symmetric nature of the stability matrix allows for a pair of complex conjugate modes to have a positive growth rate. The study of those interfacial systems has revealed characteristic differences between the nonlinear states that emerge as compared to those in related systems such as directional solidification [19] which have a regular Turing-like mode spectrum. Here the basic pattern is the traveling wave, which due to the instability extending down to k = 0, should exist up to very large wavelengths. There are thus two new features here, compared to the standard Turing case. One is the fact that the basic unstable mode is a traveling wave as opposed to a stationary spatial oscillation. It is in general quite hard to arrange parameters such as to guarantee that the basic instability is a traveling wave, while preventing a k = 0 Hopf temporally oscillating mode. Here however the marginality of the k = 0 modes makes the travellng wave instability very natural, as a k = 0 Hopf instability is impossible. Further, the standard Turing instability has a narrow band of allowed wavelengths, at least close to threshold, while here the instability extends to infinite wavelength.

It is possible to show that the traveling wave pattern arises as a supercritical bifurcation as L is increased past ${L}_{\min }=2\pi /{k}^{* }$, so that the unstable mode is an allowed wavevector. Thus the amplitude grows as $\sqrt{L-{L}_{\min }}$ as the system size is increased. For larger systems, we must turn to a numerical study of our traveling wave pattern. In the top panel of figure 2, we show an example of a traveling wave solution, corresponding to the parameter set already used above, for a periodic system of size $L=5$. The second takes the limit of infinite diffusivity for the cytosolic species c_D and c_E; for the latter case, the model is globally coupled with the value of these fields determined at all times by the integral constraints

$$\begin{eqnarray*}\begin{array}{rcl}{{Lc}}_{D} & = & {D}_{T}-\displaystyle \int {\rm{d}}{x}({c}_{d}+{c}_{{de}}),\\ {{Lc}}_{E} & = & {E}_{T}-\displaystyle \int {\rm{d}}{x}\ {c}_{{de}}\end{array}\end{eqnarray*}$$

For this size system, which is not much larger than the minimal size for the instability, ${L}_{\min }=2\pi /{k}^{* }\approx 3.0$, this traveling wave pattern appears to be the unique attractor of the system, arising from generic initial conditions.

Zoom In Zoom Out Reset image size

The solution can thus be generated by running a simulation and waiting for the system to settle into this uniformly propagating state, which perforce must be linearly stable. Alternatively, we can directly solve the steady-state equations in the moving frame of reference, so that the fields are only functions of x − vt and the equations become ordinary differential equations, by an iterative scheme acting upon the field values at collocation points. Since there are 4 second-order equations, we have to impose eight conditions. Six of these are the continuity of the four fields and two of the derivative fields across x = L. Two are the global constraints on the D_T and E_T. Because of translation invariance we can arbitrarily choose one of the fields to have a known value at say $x\,-\,{vt}=0$ and reduce the number of unknowns by one. Then the number of equations to be solved is one greater than the field unknowns, necessitating the use of the velocity as the final unknown. One can check that we get the same results from both of these methods. The second approach is specifically convenient if one has a solution for some parameter set and wishes to find a solution at a nearby one, as in that case there is a very good initial guess with which to start the iteration. We can see from these graphs that there is nothing singular about the infinite-D_c global limit at least as far as this type of solution is concerned.

We see that while the cytosol concentrations are relatively featureless (exactly so in the globally coupled, ${D}_{c}\to \infty$, limit), the membranal fields each have a single peak, located close to each other. As we take L larger, as in the bottom plot of figure 2, this structure is maintained, with the peaks having roughly the same width, and so occupying a smaller fraction of the system. The system settles into a scaling form of the solution in which the pattern consists of two parts. There is an inner region, which for the bottom frame of figure 2 lies at around the origin which gets thinner (in rescaled coordinates) as L increases. The rest of the box has an 'outer' solution which scales linearly with L. The velocity of the solution scales as L once the system is in the scaling regime, which for our parameters occurs for $L\gtrsim 15$. It is interesting to note that the large cytosol diffusion is much less successful in eliminating the variation in the cytosol fields in the larger L system.

We can understand this solution by looking separately at the two aforementioned regions, in the globally coupled limit, keeping c_D and c_E (as opposed to D_T and E_T; in the globally coupled limit, this is just a different path in parameter space, which happens to exhibit less L dependence) fixed as we increase L. We assume a dependence only on $z=(x\,-\,{vt})/L$ and hence the time derivatives become $-(v/L)\tfrac{{\rm{d}}}{{\rm{d}}{z}}$. In the outer region, the slow diffusion is irrelevant and the only spatial derivative is the velocity term. Thus, having $v\sim L$ (as can be verified from the numerical solutions) immediately allows the outer solution to have spatial decays away from the peak which become L independent in the rescaled coordinate. In the inner region, we have in general three terms that are important; the velocity term, the diffusion term and the cubic term that occurs (with opposite sign) in both the c_d and c_de equations. In fact, if we add the two equations and assume equal diffusivities, we get that ${c}_{d}+{c}_{{de}}$ does not have any driving term and one can easily check from the numerical solution that it is approximately constant on the inner scale, which is the original L independent length scale. Let us denote by A the constant value of this sum at the location of the inner zone. The numerics indicates that this quantity also scales linearly with L. If we rescale lengths by $z=\tilde{z}/L$, velocity by $v=\tilde{v}L$, and membranal concentrations by L as well, ($A=\tilde{A}L$, etc) we get the equations

$$\begin{eqnarray*}\begin{array}{rcl}0 & = & \tilde{v}\displaystyle \frac{{\rm{d}}{\tilde{c}}_{d}}{{\rm{d}}\tilde{z}}+{D}_{{\rm{m}}}\displaystyle \frac{{{\rm{d}}}^{2}{\tilde{c}}_{d}}{{\rm{d}}{\tilde{z}}^{2}}-{\lambda }_{{eE}}{c}_{E}{\tilde{c}}_{d}{(\tilde{A}-{\tilde{c}}_{d})}^{2},\\ 0 & = & \tilde{v}\displaystyle \frac{{\rm{d}}{\tilde{c}}_{{de}}}{{\rm{d}}\tilde{z}}+{D}_{{\rm{m}}}\displaystyle \frac{{{\rm{d}}}^{2}{\tilde{c}}_{d}}{{\rm{d}}{\tilde{z}}^{2}}+{\lambda }_{{eE}}{c}_{E}{\tilde{c}}_{{de}}^{2}(\tilde{A}-{\tilde{c}}_{{de}})\mathrm{.}\end{array}\end{eqnarray*}$$

This pair of conjugate equations are familiar from the literature on pattern formation, in the context of front solutions in bistable systems such as the Ginzburg–Landau equation. We can define a 'potential' function $U({\tilde{c}}_{d})$ to recast the equation for ${\tilde{c}}_{d}$, e.g. as that of a 'particle' moving with damping $\tilde{v}$ in a potential well,

$$\begin{eqnarray*}&&-\displaystyle \frac{{\rm{d}}{U}}{{\rm{d}}\tilde{z}}=\tilde{v}\displaystyle \frac{{\rm{d}}{\tilde{c}}_{d}}{{\rm{d}}\tilde{z}}+{D}_{{\rm{m}}}\displaystyle \frac{{{\rm{d}}}^{2}{\tilde{c}}_{d}}{{\rm{d}}{\tilde{z}}^{2}},\end{eqnarray*}$$

where the potential is obviously

$$\begin{eqnarray*}&&U({c}_{d})=-{\lambda }_{{eE}}{c}_{E}({\tilde{A}}^{2}{\tilde{c}}_{d}^{2}/2-2\tilde{A}{\tilde{c}}_{d}^{3}/3+{\tilde{c}}_{d}^{4}/4).\end{eqnarray*}$$

This potential has a maximum at c_d = 0 and, interestingly, a point of inflection at ${c}_{d}=A$. Since there is no minimum, for small dampling a particle starting at c_d = 0 at $\tilde{z}=-\infty$ will escape to infinite c_d. Above a critical damping, the particle approaches the point of inflection as a power-law in time, since the force is quadratic in the distance to this point. The critical damping is equal to ${\tilde{v}}_{{\rm{c}}}=\tilde{A}\ \sqrt{2{\lambda }_{{eE}}{D}_{{\rm{m}}}}$, for which the approach to the point of inflection is exponential in time, as can be seen be directly substituting in the ansatz ${c}_{d}=\tfrac{A}{2}(1+\tanh (\tilde{z}/w))$ for some width w. Thus, an inner solution exists for all velocities greater that ${\tilde{v}}_{{\rm{c}}}$. If we call this solution $\psi (\tilde{z};\tilde{v},A)$, we have the final forms ${\tilde{c}}_{d}=\psi (\tilde{z}-{\tilde{z}}_{{\rm{inner}}};\tilde{v},\tilde{A})$ and ${c}_{{de}}=\tilde{A}-{c}_{d}$. There are two unknowns, namely $\tilde{A}$ and the scaled velocity $\tilde{v}$.

The construction is then completed by integrating the L-independent outer equations in the variable z, where as mentioned above diffusion is ignored. One can then integrate the coupled first-order outer equations starting immediately past z_inner with the initial conditions ${\tilde{c}}_{d}=\tilde{A}$, ${\tilde{c}}_{{de}}=0$ and demand that the solution at $\tilde{z}=1$ returns back to the inner solution left asymptote ${\tilde{c}}_{d}=0$, ${\tilde{c}}_{{de}}=A$. These conditions determine the two unknowns. The fact that the system supports a scale-invariant nonlinear traveling wave is, we believe, traceable to the nature of the original stability; the system can use its ability to self-amplify at any non-zero k to form this solution.

While this 'single pulse' wave appears to be the unique stable steady-state solution for relatively small L, as L increases this ceases to be the case. One way to see this is to start with an initial condition which is composed of two $L/2$ pulses. For $L\lesssim 25$, for our 'standard' parameters, the two peak solution develops an instability and eventually reaches the single pulse solution appropriate to a system size of L. This process is demonstrated in figure 3 for L = 20. This is analogous to what was established long ago for a periodic array of Saffman–Taylor fingers [20]. But, here, the full story is complicated and in general we find three regions for the asymptotic state arising from these conditions. For $L\lt {L}_{c1}({D}_{{\rm{m}}})$, we obtain full coarsening as above. At larger L the two pulse pattern is stable, coexisting with the stable one pulse wave. The onset of the instability as L decreases below L_c1 is a (apparently subcritical) Hopf bifurcation. At even larger L, there is a supercritical pitchfork bifurcation to a non-symmetric two pulse state. Analogously, one can start with a three pulse initial condition in a $3L$ box and find regions of stable non-symmetric solutions at large enough L. In addition, the larger range of unstable wavevectors at increasing L appears to result in a shrinking of the basin of attraction of these steady-state solutions; this remains to be quantitatively analyzed. It should also be noted that, for finite solute diffusion, the single pulse solution ceases to exist beyond some maximal L. When this occurs, the only traveling wave state will have multiple peaks.

Zoom In Zoom Out Reset image size

Because of the stability of the single pulse solution, the system will stay in this state as the box size is adiabatically increased. To see this, we assume that the system is regulated so as to maintain a fixed overall average concentration of MinD and MinE and insert more of these proteins uniformly into the cytoplasm as we expand the cell. Following [21], the system equations then read

$$\begin{eqnarray}\begin{array}{rcl}\tfrac{\partial {c}_{D}}{\partial t} & = & \tfrac{\dot{L}(t)}{L(t)}({D}_{T}-{c}_{D})+{R}_{{de}\to D+E}-{R}_{D\to d}+\tfrac{{D}_{D}}{{L}^{2}(t)}\tfrac{{\partial }^{2}{c}_{D}}{\partial {y}^{2}},\\ \tfrac{\partial {c}_{E}}{\partial t} & = & \tfrac{\dot{L}(t)}{L(t)}({E}_{T}-{c}_{E})+{R}_{{de}\to D+E}-{R}_{d+E\to {de}}+\tfrac{{D}_{E}}{{L}^{2}(t)}\tfrac{{\partial }^{2}{c}_{E}}{\partial {y}^{2}},\\ \tfrac{\partial {c}_{d}}{\partial t} & = & -\tfrac{\dot{L}(t)}{L(t)}{c}_{d}+{R}_{D\to d}-{R}_{d+E\to {de}}+\tfrac{{D}_{d}}{{L}^{2}(t)}\tfrac{{\partial }^{2}{c}_{d}}{\partial {y}^{2}},\\ \tfrac{\partial {c}_{{de}}}{\partial t} & = & -\tfrac{\dot{L}(t)}{L(t)}{c}_{{de}}+{R}_{d+E\to {de}}-{R}_{{de}\to D+E}+\tfrac{{D}_{{de}}}{{L}^{2}(t)}\tfrac{{\partial }^{2}{c}_{{de}}}{\partial {y}^{2}},\end{array}\end{eqnarray} \tag{ 5 }$$

where $L(t)=L(0){{\rm{e}}}^{\gamma t}$ is the time-dependent length of the system and $y\equiv x/L(t)$ is the scaled spatial coordinate. Figure 4 shows clearly that the one pulse wave will maintain its global topology for a very large range of scales.

Zoom In Zoom Out Reset image size

Of course, the actual Min system does not live in a periodic domain; even if one adopts the simplification of ignoring the actual compartment structure of the cell into membrane and cytosol in favor of a bi-continuous approach (as is done here), one should obviously use zero flux conditions at the cell edges. So, it is useful to ask about the from of the nonlinear traveling wave state to the dynamics in a fixed box. In the top panel of figure 5, we present snapshots of simulations for small cells, showing clearly a 'sloshing' wave pattern with sharply decreased amplitude at the cell center. Most importantly this topology is not changed as the cell expands, even as the pattern becomes more like a traveling wave bouncing back and forth (see the bottom panel of figure 5). The time-average concentrations maintain a single node at the center even as the cell doubles; this is necessary for the functional role of the Min system in defining the precise midpoint of the cell [12, 22]. The self-adjustment property of the system allows this to take place without any fine-tuning of system parameters.

Zoom In Zoom Out Reset image size

It is critical to realize that very few of our findings should have anything to do with the detailed assumptions of the model. For example, if one uses more recent and presumably more realistic models for Min dynamics proposed in [10, 23], the existence of two conservation laws will again guarantee that the wave instability will extend down to k = 0 and therefore we can predict the existence of self-adjusting traveling wave states. This of course needs to be investigated in detail. A more uncertain situation holds for a recently studied case of a ${k}_{\min }=0$ wave instability arising during the frictional sliding of one surface above a second [24]. Here the fact that the base state with uniform sliding is explicitly not reflection symmetric and hence there need not be modes at both $+k$ and $-k$ at the same complex value of Ω; in other words, there is a preferred direction of wave propagation and this one unstable wave can be connected to just one symmetry mode as $q\to 0$. The extent to which this difference matters for the nonlinear state remains to be studied. We also expect the general features explored here will persist in a full two-dimensional model where the 2-d nature of the cytosol is treated explicitly [17].

Acknowledgments