Nonlinear self-adapting wave patterns

We propose a new type of traveling wave pattern, one that can adapt to the size of physical system in which it is embedded. Such a system arises when the initial state has an instability that extends down to zero wavevector, connecting at that point to two symmetry modes of the underlying dynamical system. The Min system of proteins in E. coli is such as system with the symmetry emerging from the global conservation of two proteins, MinD and MinE. For this and related systems, traveling waves can adiabatically deform as the system is increased in size without the increase in node number that would be expected for an oscillatory version of a Turing instability containing an allowed wavenumber band with a finite minimum.

We propose a new type of traveling wave pattern, one that can adapt to the size of physical system in which it is embedded. Such a system arises when the initial state has an instability that extends down to zero wavevector, connecting at that point to two symmetry modes of the underlying dynamical system. The Min system of proteins in E. coli is such as system with the symmetry emerging from the global conservation of two proteins, MinD and MinE. For this and related systems, traveling waves can adiabatically deform as the system is increased in size without the increase in node number that would be expected for an oscillatory version of a Turing instability containing an allowed wavenumber band with a finite minimum.
PACS numbers: 82.40.Ck,87.17.Ee,87. 16.dj 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 finite wave vector k and 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][6][7][8][9][10] To date, this instability has been viewed as an oscillatory analog of the familiar Turing instability; this means that the instability occurs at fixed q = 0, and above threshold the unstable band stretches from 0 < q min < q < q max . Here we show that there exists a new possibility, namely that for some systems with the correct symmetry, q 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 [11,12] where the static bifurcation extends to q 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 Ref. [9,13] 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 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 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 with analogous expressions for the other concentrations. This leads in a standard manner to a 4 × 4 homogeneous linear system which determines the allowed eigenvalues. The results of such a calculation for the most unstable mode are presented in Fig. 1. We see that for k < k * , there is a complex unstable mode, whose (complex) growth rate goes to 0 in the limit k → 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,3,4} and equations {2,4} 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 Ω = 0 is a doubly degenerate eigenvalue for k = 0, with left eigenvectorsφ (1) = (1, 0, 1, 1) andφ (2) = (0, 1, 0, 1) . We can then evaluate the effect of non-zero k to leading order by simply computing the projection of the diffusion terms along the diagonal, namely where the φ (j) are the corresponding right eigenvectors which satisfy the orthonomality conditionφ (i) φ (j) = δ ij . For the parameter set used in [9], the steady-state is c 0 D = 71.77, c 0 E = 76.39, c 0 d = 604.62, c 0 de = 323.61 and this gives rise to for pair-wise equal cytosol (D c ) and membrane (D m ) diffusivities. This matrix always has a pair of complex eigenvalues, as This complex pair is unstable as long as (4) which is easily satisfied by the biophysical parameters, as we saw in Fig. 1. In general, the initial rise of ReΩ 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 Ω 0 − D m k 2 for complex Ω 0 with a positive real part.
This stability structure presents a new twist on what happens in pattern forming systems such as viscous fingering and dendritic crystal growth [11,12]. There, translation invariance of the base system guarantees a single zero k = 0 eigenvalue which gives rise to a realmode instability for 0 < k < 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 [14]. 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 [15] 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 at very large wavelengths. We now turn to a study of this pattern.
In the top panel of Fig 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 The solution for the parameters of Ref. [9], for which the stability analysis is shown in Fig. 1 (blue curve). The length of the periodic system is L = 5. Middle: The solution for the limit Dc = ∞, with all other parameters as above. The membranal field profiles are basically unchanged from the above graph. Bottom: The solution for the same parameters as in the top panel, for the larger system L = 20. The peak in the membranal fields is roughly the same width as in the top panel, so that L-scaled coordinates, it appears much sharper. Away from the peak the solutal fields appear similar to the top panel, indicating that these features scale linearly with L.
The variation of the solutal fields is much increased over that of the top panel.
model is globally coupled with the value of these fields determined at all times by the integral constraints where L is the size of the periodic domain. For this size system, which is not much larger than the minimal size for the instability, L min = 2π/k * ≈ 3.0, this traveling wave pattern appears to be the unique attractor of the system, arising from generic initial conditions. 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 equation in the moving frame of reference 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 in 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 D c → ∞ 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 Fig. 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 Fig 2c lies at around x/L = 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 (data not shown), which for our parameters occurs for L 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 D and E fixed as we increase L. We assume a dependence only on z = x − vt and hence the time derivatives become −v d dz . In the outer region, the slow diffusion is irrelevant and the only spatial derivative is the velocity term. Thus, having v ∼ L 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 a smooth function on the inner scale. Let us denote by A the constant value of this sum at the location of the inner zone. If we rescale lengths by z =z/L, velocity by v =ṽL, we get the equations This pair of conjugate equations are familiar from the literature on pattern formation. We can define a "potential" function U (c d ) to recast the equation for c d , e.g. as that of a "particle" moving a well where the potential is obviously . The inner solution is then a particle which rolls as z goes from −∞ to +∞ from the U maximum at c d = 0 to the U point of inflection at c d = A. From this analogy, it is obvious that solutions of this form exist for all velocities above a critical velocity which can be found to equal A √ 2λ eE Dm L by directly substituting in the ansatz c d = A 2 (1 + tanh(z/w)). Note that for velocities above critical, there is actually a power-law decay to the fixed point at A rather than the exponential decay obtained for the minimal velocity. If we call this solution ψ(z; v, A), we have the final forms c d = ψ(z − z inner ;ṽ, A) and c de = A−c d . There are two unknowns, namely A and the scaled velocityṽ.
The construction is then completed by integrating the outer equations, where diffusion is ignored and lengths are now scaled as (z − z inner ) =zL; this choice cancels the L dependence in the velocity term and gives us an equation where all terms are O(1) and so L has completely disappeared from the problem. One can then integrate the coupled first-order outer equations starting immediately past z inner with the initial conditions c d = A, c de = 0 and demand that the solution atz = 1 The only possible difficulty is that the determined velocityṽ could fall below the minimum velocity of the inner solution; since that latter scales as 1/L this will always occur as L is decreased and in fact defines the lower L limit for the existence of this self-adjusting wave. Waves do continue to exist below this L, but they have a more complex scaling indicated for example by the amplitude A changing with L. 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 q 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 25, for our "standard" parameters, the two peak solution develops an instability and eventually reaches the single peak solution appropriate to a system size of L. This process is demonstrated in Fig 3 for L = 20. This is analogous to what was established long ago for a periodic array of Saffman-Taylor fingers [16]. But, here, the full story is complicated and in general we find three regions for the asymptotic state arising from these conditions. For L < L c1 (D m ), we obtain full coarsening as above; at larger L the two pulse pattern is stable, coexisting with the stable one pulse wave. 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.
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 Ref. [17]. the system equations then read where L(t) = L(0)e γt is the time-dependent length of the system and y ≡ x/L(t) is the scaled spatial coordinate. Fig. 4 shows clearly that the one pulse wave will maintain its global topology for a very large range of scales. 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 Fig. 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 Fig. 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 [10,18]. The self-adjustment property of the system allows this to take place without any fine-tuning of system parameters.
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 refs. [8,19], the existence of two conservation laws will again guarantee that the wave instability will extend down to q = 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 q min = 0 wave instability arising during the frictional sliding of one surface above a second [20]. Here the fact that the base state with uniform sliding is explicitly not reflection symmetric and hence there need not be modes at both +q and −q 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 → 0. The extent to which this difference matters for the non-linear state remains to be studied. with reflecting boundary conditions. Bottom: The adiabatic adjustment of the above "sloshing" pattern in an exponentially growing system, with a doubling time of 30 mins. and an initial size L = 1.5 The system quickly settles into the sloshing pattern, which for larger systems resembles the onepulse traveling wave solution when the pulse is away from the end-walls. The other parameters are the same as in the top panel of Fig. 2.