1. Introduction

One of the most striking features of living systems is the appearance of seemingly complex, sophisticated patterns that often evolve with time and/or the development of the organism. Such patterns are most obvious to humans at the macroscopic level, for example, animal coat patterns or patterns in flowers and seeds. These patterns follow sophisticated mathematical formulae, symmetries and fractal dimensions that have puzzled scientists for millennia. Underlying these patterns are self-organising molecular systems where local interactions create global structure.

Inspired by biological patterns, Turing developed a mathematical model that would describe spontaneous pattern formation based on reaction-diffusion equations [1]. Within this work, he showed that the addition of diffusion terms to non-linear reactions can result in a destabilisation of the spatially homogeneous steady state giving rise to spatial patterning. Similarly, Gierer and Meinhardt discovered sets of reaction-diffusion equations based on short-range activation and long-range inhibition that produced spatiotemporal patterning [2]. At the organismal level, these approaches have been used to produce theoretical models for developmental biology [3].

Patterning is also a characteristic feature of single cells. Many proteins and cellular structures reside at particular cellular sites that are related to the geometry of the cell. Examples include division septa, flagella, pili, and stalks in bacteria [4]. Since the 1990s, it has become clear that for many bacterial proteins, accurate localization is essential for correct function. Examples include: the clustering of chemotactic receptors at the poles of cells [5,6]; the compartment-specific production and assembly of sporulation proteins governing spore morphogenesis [7,8,9]; and the mechanism of cell division [10,11].

The accurate localisation of proteins and subcellular structures to specific sites in a cell predicates the existence of physical mechanisms capable of determining the positions of sites within cells before locating the appropriate proteins. In bacteria, the most studied example of a subcellular localization mechanism is the positioning of the FtsZ ring for cell division in Escherichia coli [12]. The Min protein system is involved in establishing the timing and position of the FtsZ contractile ring that ultimately separates the cell into two daughter cells [13]. Prior to cell division, membrane-bound Min proteins form an oscillating pattern where the proteins are concentrated at one pole of the cell and then another, leaving a bare zone or node at the cell equator where the FtsZ division ring will form [13].

In its simplest form, the Min system consists of two proteins, MinD and MinE, a lipid bilayer (membrane) and ATP [14]. MinD is an ATPase that exists in two forms: a soluble form and a membrane-bound form. When bound to ATP, MinD binds to the lipid bilayer. MinE is a soluble protein that can bind to membrane-bound MinD.ATP and accelerate the hydrolysis of the MinD-bound nucleotide. The membrane-bound MinE-MinD.ADP complex then dissociates, returning all components into solution. Soluble MinD rapidly rebinds ATP and returns to the membrane-bound state.

In vitro, this system has been shown to spontaneously produce travelling spiral wave patterns of membrane-bound MinD and MinD.MinE on a planar lipid bilayer [14]. The moving spiral waves resemble those seen in Belousov-Zhabotinsky reactions [15]. The observation of spiral wave patterns produced by purified proteins in an artificial system indicates that in vivo these four components (MinD, MinE, ATP and lipid bilayer) are all that is required for pattern formation. Clearly other components are required to locate the FtsZ ring at the node in oscillating Min pattern on the cell membrane, but these are not likely to be required for establishing the Min pattern.

Several studies have modelled the Min protein system using reaction-diffusion equations [16,17,18,19]. A recent partial differential equation model based on physical protein interactions accurately reproduces the spatiotemporal patterns that correspond to the observed Min oscillations in cells [20]. In these models, pattern formation arises from an instability in the reaction-diffusion equations. These instabilities are essentially “Turing patterns” [1], where the non-linear reactions couple with diffusion to destabilise an otherwise homogeneous steady state.

Systems that lead to pattern formation within cells are often poorly characterized or their components completely unidentified. As a result, it is impossible to construct complete mathematical models that give rise to the observed patterns. Faced with this lack of information, it is still possible to derive heuristic, top down, models that can predict some of the observed features of the pattern forming system.

There is a link between the morphology of a cell and the formation of patterns within the cell [21]. By assuming the existence of an underlying reaction-diffusion system, we aim to develop a set of mathematical tools that elucidate such a link. Whilst the patterning of the Min system in E. coli is well understood, such tools may be useful in detecting and characterising other intracellular patterning systems, such as those that must underlie cell division in other organisms. To verify our heuristic model, we will use the Min system as an example. This allows us to compare the heuristic model to both a realistic reaction-diffusion mathematical model as well as direct observations of Min patterning in E. coli.

In this paper, we show that the general features of the spatiotemporal patterning of the Min protein system are well modelled using concepts that originate from a general approach to linear stability analysis of reaction-diffusion equations. As part of this linear stability analysis a dispersion relation is calculated that relates the stability of the homogeneous steady state to the eigenspectrum of the diffusion operator, or spatial frequencies [22]. The linear stability analysis of a reaction-diffusion model of the Min system identifies the presence of the Turing instability, but fails to quantitatively predict other geometry-dependent features of the dynamics such as the presence of oscillations and the dominant wavenumber (or spatial frequency) of the pattern.

We construct a heuristic model by assuming that the principles that lead to the dispersion relation hold in a more general sense. The effects of reaction dynamics and cell shape can then be separated, with the link between them expressed through a heuristic dispersion relation. This heuristic dispersion relation can be parameterised by observing patterns on a range of cell geometries or sizes. Once the heuristic dispersion relation is parameterised it should be able to predict the observed features of the pattern on any geometry without any knowledge of the true underlying reactions. The model will predict that the patterns should approximately correspond to harmonic modes or eigenfunctions of the Laplace-Beltrami operator over the cell shape.

To test the validity of our heuristic approach we show that the heuristic model can approximate the spatial distribution of the Min proteins from the geometry of the cell without any knowledge of the underlying dynamics. To do this, we derive a heuristic dispersion relation from a set of full solutions of our Min model [20] calculated on a growing, filamentous cell that does not divide. As the cell grows, the spatial patterns of Min protein distributions follow a sequential transition through the harmonic modes and this data is used to parametrise the heuristic dispersion relation. We use this heuristic dispersion relation to predict the patterning of the Min system on an arbitrary cell shape and show that the resulting MinD protein distributions compare favourably with complete solutions of our full theoretical model. Subsequently, we show that the heuristic model predicts patterns that are consistent with published experimental data for E. coli cells that have been deformed to arbitrary shapes via nanofabrication [23,24].

2. Materials and Method

2.1. Model overview

The model of the Min system used in this work is adapted from a model previously proposed by Walsh et al. [20]. This model contains six species: MinD monomer in solution (D), MinE homodimer in solution (E₂), membrane-bound MinE homodimer (e₂), membrane-bound MinD monomer (d) and dimer (d₂), and a membrane bound MinD/E heterotetramer complex (d₂e₂).

This molecular model is based on the experimentally identified interactions of the Min proteins. In particular, crystal structures have shown that MinD monomers (d) are able to form dimers (d₂) in an ATP-dependent manner [25]. Furthermore, FRET and yeast two hybrid experiments have shown that this occurs in a two-step process. MinD monomers in solution (D) first bind to the membrane as a monomer (d) before forming a dimer while bound to the membrane (d₂) [26,27]. MinE is a stable homodimer which undergoes a large-scale structural change between its inactive cytosolic and active membrane-bound forms [28,29,30,31].NMR experiments have shown that MinE heavily favors its inactive cytosolic form (E₂) [30] unless it is stabilized in its active conformation by binding to a membrane-bound MinD dimer to form a MinD/E complex (d₂e₂) [28]. While in this complex, MinE allosterically activates MinD to hydrolyze ATP [32] causing the complex to dissociate with MinD being released into the cytosol as monomers (D) while MinE returns to being an unstable membrane-bound homodimer (e₂).

In order for the molecular model of the Min system to function, it was hypothesized that membrane-bound MinD monomer (d) can be removed from the membrane by interacting with membrane-bound MinE (e₂). The inclusion of this interaction allows the model to recreate many of the characteristics of Min patterning seen in vivo [20].

These experimentally determined properties of the Min proteins are incorporated into the model of the Min system by the following reactions: MinD monomer in solution (D) can bind to the membrane to give rise to a membrane-bound MinD monomer (d). The membrane-bound MinD monomer (d) can either form a dimer (d₂), or be released from the membrane where the latter is stimulated by the membrane-bound MinE homodimer (e₂). The MinD dimer (d₂) can bind to the membrane-bound MinE (e₂) to form a MinD/E heterotetramer complex (d₂e₂). Following ATP hydrolysis, this complex dissociates resulting in MinD being released into the cytosol as monomers (D) while MinE returns to being an unstable membrane-bound homodimer (e₂). MinE dimer in solution (E₂) is able to bind to the membrane to give rise to the membrane-bound homodimer (e₂) which, in turn, can be released back into the cytosol. These simplified reactions are shown schematically in Figure 1.

To simplify the model of the Min system, we assume that the cytosolic states are homogenous perpendicular to the membrane. As a result of this approximation, the value of the cytoplasmic states (MinD in solution (D) and MinE in solution (E₂)) only needs to be measured on the membrane, reducing the dimensionality of the model by one. Further, all changes to the cytoplasmic states through binding and release from the membrane will be proportional to the surface-to-volume ratio of the shape (S). Solving the model on a one-dimensional line with a surface to volume ratio of $4\mu {m^{ - 1}}$ (the value used for all simulations in this work) is approximately equivalent to solving the original equations on a cylinder with a radius of $0.5\mu m$.

The resulting set of PDEs in this work are,

$\begin{array}{l} {\partial _t}D = S\left( {{\omega _{edf}}{e_2} \cdot d + 2{\omega _{hydr}}{d_2}{e_2} - {\omega _{db}}D} \right) + {D_D}{\nabla ^2}D, \\ {\partial _t}{E_2} = S\left( {{\omega _{er}}{e_2} - {\omega _{eb}}{E_2}} \right) + {D_E}{\nabla ^2}{E_2}, \\ {\partial _t}d = {\omega _{db}}D - 2{\omega _{dim}}{d^2} - {\omega _{edf}}{e_2} \cdot d + {D_m}{\nabla ^2}d, \\ {\partial _t}{d_2} = {\omega _{dim}}{d^2} - {\omega _{edf}}{e_2} \cdot {d_2} + \frac{1}{2}{D_m}{\nabla ^2}{d_2}, \\ {\partial _t}{e_2} = - {\omega _{er}}{e_2} + {\omega _{eb}}{E_2} - {\omega _{edf}}{e_2} \cdot {d_2} + {\omega _{hydr}}{d_2}{e_2} + \frac{1}{2}{D_m}{\nabla ^2}{e_2}, \\ {\partial _t}{d_2}{e_2} = {\omega _{edf}}{e_2} \cdot {d_2} - {\omega _{hydr}}{d_2}{e_2} + \frac{1}{4}{D_m}{\nabla ^2}{d_2}{e_2}. \end{array}$

The rates and diffusion coefficients adapted from Walsh et al. [20] can be found in Table 1.

2.2. Construction of the dispersion relation

The PDE describing the model of the Min system can be rewritten as ${\partial _t}{\rm{c = R + D}}{\nabla ^2}{\rm{c}}$ using three vectors: the species vector containing the concentrations of each species in order (${\rm{c}} = \left\{ {D, {E_2}, d, {d_2}, {e_2}, {d_2}{e_2}} \right\}$), the diffusion vector which contains the diffusion constants for the species in the same order (${\rm{D = }}\left\{ {{D_D}, {D_E}, {D_m}, {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}{D_m}{\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}{D_m}, {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 4$}}{D_m}} \right\}$), and the reaction vector, containing the reactions generating each species in the same order:

${\rm{R = }}\left( {\begin{array}{*{20}{c}} {S\left( {{\omega _{edf}}{e_2} \cdot d + 2{\omega _{hydr}}{d_2}{e_2} - {\omega _{db}}D} \right)}\\ {S\left( {{\omega _{er}}{e_2} - {\omega _{eb}}{E_2}} \right)}\\ {{\omega _{db}}D - 2{\omega _{dim}}{d^2} - {\omega _{edf}}{e_2} \cdot d}\\ {{\omega _{dim}}{d^2} - {\omega _{edf}}{e_2} \cdot {d_2}}\\ { - {\omega _{er}}{e_2} + {\omega _{eb}}{E_2} - {\omega _{edf}}{e_2} \cdot {d_2} + {\omega _{hydr}}{d_2}{e_2}}\\ {{\omega _{edf}}{e_2} \cdot {d_2} - {\omega _{hydr}}{d_2}{e_2}} \end{array}} \right)$

The homogeneous steady state is found by solving R = 0 with the additional restrictions that the total concentration of MinD is set (${\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 4$}}D + d + 2{d_2} + 2{d_2}{e_2} = 347.25\mu {m^{ - 2}}$) and the total concentration of MinE is set (${\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 4$}}{E_2} + {e_2} + {d_2}{e_2} = 121.5\mu {m^{ - 2}}$). The resulting homogenous steady state species vector ${\rm{c = }}\left\{ {30.7\mu {m^{ - 3}}, 22.2\mu {m^{ - 3}}, 83.4\mu {m^{ - 2}}, 12.2\mu {m^{ - 2}}, 0.05\mu {m^{ - 2}}, 115.9\mu {m^{ - 2}}} \right\}$ is which gives the concentrations of each species.

A linear approximation for the behaviour of the system around the steady state can be obtained by evaluating the Jacobian matrix of the reaction vector that is given by ${{\rm{J}}_{ij}} = \frac{{\partial {{\rm{R}}_i}}}{{\partial {{\rm{c}}_j}}}$, which for this model is:

${\rm{J}} = \left( {\begin{array}{*{20}{c}} { - S{\omega _{db}}}&0&{S{\omega _{edf}}{e_2}}&0&{S{\omega _{edf}}d}&{2S{\omega _{hydr}}}\\ 0&{ - S{\omega _{eb}}}&0&0&{S{\omega _{er}}}&0\\ {{\omega _{db}}}&0&{ - 4{\omega _{{\rm{dim}}}}d - {\omega _{edf}}{e_2}}&0&{ - {\omega _{edf}}d}&0\\ 0&0&{2{\omega _{{\rm{dim}}}}d}&{ - {\omega _{edf}}{e_2}}&{ - {\omega _{edf}}{d_2}}&0\\ 0&{{\omega _{eb}}}&0&{ - {\omega _{edf}}{e_2}}&{ - {\omega _{er}} - {\omega _{edf}}{d_2}}&{{\omega _{hydr}}}\\ 0&0&0&{{\omega _{edf}}{e_2}}&{{\omega _{edf}}{d_2}}&{ - {\omega _{hydr}}} \end{array}} \right)$

To construct the dispersion relation we evaluate the Jacobian matrix at the homogeneous steady state and numerically determine the eigenvalues of the matrix ${\rm{M}} = {\rm{J}} - {k^2}{\rm{ID}}$ (where I is the identity matrix). The eigenvalue with the largest real component is plotted against the parameter k. This is seen in Figure 2.

2.3. Numerical methods

2.3.1. Calculation of the Dispersion Relation

The numerical calculations for the construction of the dispersion relation were carried out using Mathematica 10.2. The homogeneous steady state (solution to R = 0) was found by numerically solving the system of ODEs

$\frac{{d{\rm{c}}}}{{dt}} = {\rm{R}}$

for 1000 seconds using the initial conditions ${\rm{c = }}\left\{ {1389\mu {m^{ - 3}}, 486\mu {m^{ - 3}}, 0, 0, 0, 0} \right\}$ (i.e. starting with all protein in the cytoplasm). The ODE was solved using the NDSolve function. The run time of 1000 seconds was sufficient to show no further variations, and the solutions were taken as numerical approximations for the steady state values. This method ensured that the total concentration of proteins remained fixed and that the solution found was stable.

The dispersion curve was calculated using the matrix ${\rm{M}} = {\rm{J}} - {k^2}{\rm{ID}}$ by setting k to values between zero and 4.5 in steps of 0.01. At each step, the eigenvalues of the matrix were found using the Eigenvalues function. The eigenvalue with the largest real component at each step was plotted versus the value of k to give the dispersion curve.

2.3.2. One dimensional simulations

The one spatial dimension nonlinear PDE model was solved using Mathematica by employing the method of lines. A structured point centre mesh consisting of 100 points was manually generated to form a set of 600 coupled ODE equations which were solved using NDSolve. Initial conditions were taken as the homogeneous steady state concentrations multiplied by a uniformly distributed random value between 0.95 and 1.05 at each node point.

Mode of patterning was calculated by first running the simulation for 1000 seconds. The largest discrete cosine transform component of the summed MinD concentration was found by applying the FourierDCT function at each second over a 250 second period. The most common maximum component over the 250 seconds was taken to be the dominant mode.

The random noise applied to the initial value of can cause the system to pattern in an energetically unfavourable mode. To find the most energetically favourable and hence most probable mode, 10 repeats at each length were simulated and the value of the most common dominant mode retained. The resulting mode transition plots are shown in Figure 4B and 4C.

2.3.3. Eigenfunctions of two-dimensional shapes

The eigenfunctions and eigenvalues of the Laplace-Beltrami operator on two dimensional shapes were found using Mathematica. The shape was meshed by first creating a Polygon Object using the Polygon function on the outline of the shape. Then, the DiscretizeRegion function was applied to the object to create the mesh. The NDEigensystem function was used on the mesh of the shape, with zero-flux boundary conditions enforced, to calculate its eigenfunctions.

2.3.4. Solution of the non-linear PDE model in two dimensions

The solutions of the non-linear PDEs on two-dimensional shapes were computed using FlexPDE (version 6.32, PDE Solutions Inc.). Parameters, concentrations and initial conditions were as per the one-dimensional simulations. Simulations were run for 1000 seconds before the patterning was analysed every 10 seconds over a period of 250 seconds. Length steps of $0.25\mu m$ between simulations were taken. The dominant mode of each simulation was calculated using Mathematica by finding the eigenfunction with the absolute maximum inner product between itself and the simulation data.

2.4. Two dimensional shapes formation

The arbitrary shape used for Figure 6 was created by randomly generating 20 uniformly distributed random values between −0.4 and 0.4, which were then used as the coefficients for a five Fourier basis expansion of a parametric outline. Numbers were redrawn until a non-intersecting outline was formed.

3. Results

3.1. Theoretical framework: Coordinate free linear stability analysis

Traditionally, linear stability analysis is performed on geometries with Cartesian symmetry, be it a one-dimensional line, a rectangle or a rectangular prism. However, it is possible to take a generalized approach to show that linear stability can be used to predict patterning on any shape or volume. Such a coordinate free approach to linear stability analysis has been proposed previously [33]. A derivation with consistent notation for this paper is presented below. Readers who are less mathematically inclined can skip this section.

In general, a reaction-diffusion model can be described by the set of equations,

$\frac{{\partial {\rm{c}}\left( {{\rm{x}}, t} \right)}}{{\partial t}} = {\rm{R}}\left( {{\rm{c}}\left( {{\rm{x}}, t} \right)} \right) + {\rm{D}}\Delta {\rm{c}}\left( {{\rm{x}}, t} \right)$

(1)

where ${{\rm{c}}\left( {{\rm{x}}, t} \right)}$ is a vector of spatially and temporally dependent concentrations for each of the constituent species, ${\rm{R}}\left( {{\rm{c}}\left( {{\rm{x}}, t} \right)} \right)$ is the non-linear reaction vector that describes the interaction between species, D is the vector of diffusion constants for each species, and $\Delta$ is the linear operator that describes the diffusive motion. These equations operate in some domain, M, which in this instance describes the geometry of the cell. If the domain M is in a Euclidean plane, $M \subset {{\rm{R}}^2}$ then $\Delta = {\nabla ^2}$ is the Laplacian operator. More generally if M is a Riemannian manifold then $\Delta {\rm{c = div}}\left( {{\rm{grad}}\left( {\rm{c}} \right)} \right)$ is the Laplace-Beltrami operator. In general the reaction-diffusion equations can only be solved numerically. However, the spatiotemporal behaviour of the system can be determined by a semi-analytic approach by taking a linear expansion around a homogenous equilibrium solution.

Provided that the homogenous steady state is stable in the absence of diffusion, i.e. with D = 0, then linearizing the dynamics around this point will provide an approximation of the full dynamics including diffusion. Labelling the steady state concentrations as c₀, and writing deviations from this steady state as $\delta {\rm{c}}\left( {{\rm{x}}, t} \right) = {\rm{c}}\left( {{\rm{x}}, t} \right) - {{\rm{c}}_0}$ the linearized system is given by,

${\delta _t}\left( {\delta {\rm{c}}\left( {{\rm{x}}, t} \right)} \right) = {\rm{J}}\left( {{{\rm{c}}_0}} \right) \cdot \delta {\rm{c}}\left( {{\rm{x}}, t} \right) + {\rm{D}}{\nabla ^2}\delta {\rm{c}}\left( {{\rm{x}}, t} \right),$

(2)

where ${\rm{J}}\left( {{{\rm{c}}_0}} \right)$ is the linearization of the reactions vector calculated by taking the Jacobian of ${\rm{R}}\left( {{\rm{c}}\left( {{\rm{x}}, t} \right)} \right)$, evaluated at the homogenous steady state, so that ${\rm{R}}\left( {{{\rm{c}}_0} + \delta {\rm{c}}\left( {{\rm{x}}, t} \right)} \right) \approx {\rm{J}}\left( {{{\rm{c}}_0}} \right) \cdot \delta {\rm{c}}\left( {{\rm{x}}, t} \right)$. Deviations from the steady state can be written as a sum of separable solutions, $\delta {\rm{c}}\left( {{\rm{x}}, t} \right) = {\Sigma _k}{\rm{X}}\left( {{\rm{x}}, k} \right){\rm{T}}\left( {t, k} \right)$, so long as the set of spatial functions, ${\rm{X}}\left( {{\rm{x}}, k} \right)$, span the Hilbert space L² over the domain of the system (the cell). This becomes useful if each term in the sum satisfies the linearized equation as then each basis component becomes decoupled and can be examined independently. Substituting a single term of the separable solution into the linearized equation (Equation 2) and solving yields,

$\frac{{{\partial _t}\left( {{\rm{T}}\left( {t, k} \right)} \right) - {\rm{J}}\left( {{{\rm{c}}_0}} \right) \cdot {\rm{T}}\left( {t, k} \right)}}{{{\rm{D}} \cdot {\rm{T}}\left( {t, k} \right)}} = - {k^2} = \frac{{{\nabla ^2}{\rm{X}}\left( {{\rm{x}}, k} \right)}}{{{\rm{X}}\left( {{\rm{x}}, k} \right)}}.$

(3)

From the right hand side of this equation, the spatial functions, ${{\rm{X}}\left( {{\rm{x}}, k} \right)}$, that satisfy the linearized system are eigenfunctions of the Laplace-Beltrami operator, also known as harmonic modes. Each of the eigenfunctions will be coupled to the reactions through its eigenvalue, -k².

The time evolution of the amplitude of each of the spatial basis function over time is calculated by solving the left hand side of Equation 3, which yields,

${\rm{T}}\left( {t, k} \right) = {\rm{T}}\left( {0, k} \right) \cdot {e^{{\rm{M}}\left( k \right)t}},$

(4)

where the matrix ${\rm{M}}\left( k \right) = {\rm{J}}\left( {{{\rm{c}}_0}} \right) - {k^2}{\rm{ID}}$ is a function of the separation constant -k². The solution can be expressed as,

${\rm{T}}\left( {t, k} \right) = \sum\limits_{i = 0}^n {{e^{{\lambda _i}\left( k \right)t}}} {{\rm{P}}_i}\left( {t, k} \right){{\rm{V}}_i}\left( k \right)$

(5)

where ${{\lambda _i}\left( k \right)}$ is the i^th eigenvalue of ${\rm{M}}\left( k \right)$, ${{\rm{V}}_i}\left( k \right)$ is the i^th eigenvector of ${\rm{M}}\left( k \right)$ and the polynomial ${{\rm{P}}_i}\left( {t, k} \right)$ is dependent on the degeneracy of the i^th eigenvalue of ${\rm{M}}\left( k \right)$ and the initial condition of ${\rm{T}}\left( {0, k} \right)$.

Due to the exponentials in equation (5), at large values of t, ${\rm{T}}\left( {t, k} \right)$ will be dominated by a single term in the sum. This term corresponds to the eigenvalue with the largest real component, which we will denote ${\lambda _{\max }}\left( k \right) = \mathop {\max }\limits_i \left\{ {{\mathop{\rm Re}\nolimits} \left( {{\lambda _i}\left( k \right)} \right)} \right\}$.

The perturbation, $\delta {\rm{c}}\left( {{\rm{x}}, t} \right) = {\Sigma _k}{\rm{X}}\left( {{\rm{x}}, k} \right){\rm{T}}\left( {t, k} \right)$, can only grow in time if one of the time dependent functions, ${\rm{T}}\left( {t, k} \right)$, grows in time. Thus, the sign of the eigenvalue ${\lambda _{\max }}\left( k \right)$ dictates the linear stability of the system. If for all k then the system is stable, as all perturbations are damped. A Turing instability is defined by ${\lambda _{\max }}\left( 0 \right) < 0$ and ${\lambda _{\max }}\left( k \right) > 0$ for some k > 0 [1,22]

3.2. Theoretical results

3.2.1. Linear stability analysis provides a method to separate reactions from geometry

When checking for a Turing instability in a system of reaction-diffusion equations, the results of a linear stability analysis are typically visualized using a dispersion relation (Figure 2). For a dispersion relation, the curve can be found by considering only the reaction terms and the diagonal matrix of diffusion coefficients. The curve encodes the stability of the system as a function of a parameter k, the wavenumber, which has the units of inverse length. When the curve is negative at and positive for some range of k values then the system of equations has a Turing instability.

The dispersion curve is independent of the geometry of the system, e.g. the shape of the cell. However, the geometry dictates the sampling of the dispersion curve. When the system is restricted to a bounded domain the k parameter may only take values that are related to the eigenvalues of the Laplace-Beltrami operator on that domain. For each permitted value of k there will be a corresponding eigenfunction, or harmonic mode, that may approximate the spatial distribution of particles in the system (MinD in our case).

For the case of a one-dimensional line, the Laplace-Beltrami operator is ${\nabla ^2}{\rm{c}}\left( {{\rm{x}}, t} \right) = \frac{{{\partial ^2}}}{{\partial {x^2}}}{\rm{c}}\left( {{\rm{x}}, t} \right)$, with eigenfunctions of the form $c\left( {x, t} \right) = a\sin \left( {kx} \right) + b\cos \left( {kx} \right)$, were a and b are arbitrary constants. Differentiating the eigenfunction twice confirms that the associated eigenvalue is -k². When restricted to a finite domain of length l, subject to zero flux boundary conditions, the eigenfunction must obey the boundary condition. This in turn restricts the allowed values of the wavenumber, k, and the zero flux boundaries give $k = \frac{{n\pi }}{l}$, where n is a non-negative integer. Each of these wavenumbers correspond to a harmonic mode where the reactant distribution is given by $\cos \left( {\frac{{n\pi }}{l}x} \right)$ where n is the spatial mode number.

The dispersion relation shown in Figure 2 was calculated by applying linear stability analysis to the reaction-diffusion equations used previously to accurately model the E. coli Min system [20] (see Methods for details). It shows a region where ${\lambda _{\max }}\left( k \right)$ is positive ($0.18\mu {m^{ - 1}} < k < 3.35\mu {m^{ - 1}}$), thus, if the cell geometry produced a wavenumber, k, in this range, the corresponding eigenfunction will be amplified, resulting in spatial patterning of Min proteins. The dispersion curve is negative for $k < 0.18\mu {m^{ - 1}}$ and $k > 3.35\mu {m^{ - 1}}$ thus eigenfunctions whose corresponding wavenumbers lie in these regions will be damped. Critically, the dispersion curve is negative at k = 0, and hence this system obeys the conditions for a Turing pattern.

The geometry and size of the cell dictates the sampling of the dispersion curve. In the example shown in Figure 2, the cell is represented by a one-dimensional object, of length $3.35\mu m$. This results in a sampling of the dispersion curve at the wavenumbers $k = 0, 0.94, 1.88, 2.81{\rm{ and }}3.75\mu {m^{ - 1}}$ (five blue points) representing the lowest five eigenfunctions. As stated previously, the dispersion curve is negative for the first eigenvalue, k = 0. The dispersion curve is positive for the next three wavenumbers $k = 0.94, 1.88, {\rm{ and }}2.81\mu {m^{ - 1}}$, thus, the modes corresponding to these values will be amplified. It can be assumed that the eigenfunction corresponding to the most positive value of the dispersion curve will be the most heavily amplified and will dominate the behaviour of the system.

For a longer cell the spacing between the sampling points will be smaller, as the k values have units of inverse length. Conversely, for shorter cell lengths, the spacing between these points will increase.

There are two key points that are illustrated by a coordinate free approach to linear stability analysis. Firstly, the dispersion curve is independent of the geometry of the domain. The geometry only dictates the points on the dispersion curve that are sampled. The resulting patterning will obviously change depending on the shape of the domain. However, the dispersion curve and the wavenumbers, k, at which the reactions support patterning will not change.

The second key point is that regardless of what reactions underlie a Turing pattern, linear stability analysis predicts that the pattern will approximate an eigenfunction of the Laplace-Beltrami operator. In geometries where the Laplace-Beltrami eigenvalues are degenerate, or large geometries where there are multiple eigenfunctions that are strongly amplified, then mode mixing may occur and the pattern may not be well approximated by a single eigenfunction. However, in small domains where the sampling wavenumbers are well spaced, the pattern is likely to be well approximated by a single eigenfunction.

As a cell grows its geometry will change, leading to a change in the permitted wavenumbers. This linear approximation would then predict that the pattern should also change as the cell grows. For a one-dimensional line, ignoring possible hysteresis effects, as the size of the domain is increased the pattern is expected to transition through each of the eigenfunctions.

3.2.2. Linear stability analysis correctly predicts the onset of patterning for the Min protein system

Patterning occurs when an eigenvalue of the Laplace-Beltrami operator over the geometry of the cell samples the region of the dispersion curve that is positive. Thus, the system will not exhibit patterning when the lowest non-zero eigenvalue lies beyond this positive region of the dispersion curve. For the example in Figure 2, this occurs when the wavenumber $k = 3.346\mu {m^{ - 1}}$, which corresponds to the lowest order non-zero eigenvalue for a cell length of $0.938\mu m$. Thus, linear stability analysis predicts that for cells shorter than this length, spatial MinD patterns will not be observed.

This prediction can be compared to the results of a full simulation of the non-linear equations describing the Min protein system. These are visualized using kymograms (Figure 3), where the concentration of MinD along the one-dimensional line (corresponding to the long axis of the cell) is plotted on the y-axis with time running along the x-axis. High concentrations of MinD are coloured in yellow while low MinD concentrations are dark blue.

As shown in Figure 3A at a length of $0.93\mu m$ no patterning occurs for the full simulation of the non-linear equations (i.e. the MinD distribution is homogeneous) while at $0.94\mu m$ a stable stationary first order mode is created. Thus, linear stability analysis correctly predicts the onset of patterning for our model of the Min protein system.

3.2.3. Linear stability analysis fails to predict the onset of Min oscillations

Oscillating patterns can be predicted from dispersion curves. They require the eigenvalues of the evolution operator, ${\lambda _{\max }}\left( k \right)$, to be complex and have a positive real component. In Figure 2, for the model of the Min protein system, the regions of the dispersion curve where ${\lambda _{\max }}\left( k \right)$ is real are shown in green, while the regions where it is complex are shown in red. Based on this, the linear stability analysis of this system predicts that the MinD patterns should be stationary for all cell lengths. However, full simulation of the non-linear equations clearly show that MinD patterns oscillate (Figure 3B) with the onset of oscillations occurring at a cell length of $2.7\mu m$, which corresponds to a wavenumber $k = 1.16\mu {m^{ - 1}}$. Thus, while the full simulation of the non-linear equations reproduces the experimentally observed MinD oscillations, the linear stability analysis of the same equations does not.

3.2.4. Linear stability analysis qualitatively predicts the transitions to higher order modes

Linear stability analysis can be used to determine when the MinD pattern will move from a first order mode to a second order mode (i.e. MinD distributions described by $\cos \left( {\frac{\pi }{l}x} \right)$ and $\cos \left( {\frac{{2\pi }}{l}x} \right)$ where x is the position on the long axis of the cell and l is the cell length). This is expected to occur when the ${\lambda _{\max }}\left( k \right)$ values for the two corresponding eigenvalues are equal, hence, equally amplified. In Figure 2, we see that for a cell length of $3.35\mu m$ both the first and second modes have the same level of amplification (i.e. same value of ${\lambda _{\max }}\left( k \right)$). If the length of the cell were slightly shorter the sampling wavenumbers would be shifted to the right. As a result, the first order mode would be more heavily amplified than the second and so linearly stability analysis would predict that the spatial distribution of MinD would approximate the first harmonic mode. The reverse is true if the length were increased, favouring the second order mode.

By simulating the full non-linear equations describing the Min system, we have shown that the system transits from the first order mode to the second order mode at a cell length of $5.5\mu m$. This is shown in Figure 4A (red curve) where the edge of the first step occurs at $5.5\mu m$. The red curve in Figure 4A shows the cell lengths at which the MinD pattern transits through successive modes as determined by full simulation of the non-linear equations describing the system. These can be compared with the cell lengths for the same transitions as predicted by linear stability analysis of the same equations (Figure 4A, blue curve). As shown by the curves, linear stability analysis consistently predicts the mode transitions to occur before they are actually observed in the full simulation.

3.2.5. Heuristic dispersion relation derived from the behavior of the Min system

Despite the failings of linear stability analysis to predict the Min oscillations and to predict the cell lengths at which the distribution of MinD transits between consecutive modes, many of the core tenets of the linear stability analysis still hold true.

The spatial patterning of MinD closely approximates the eigenmodes dictated by the cell geometry. Throughout the various full simulations of the non-linear equations (Figure 3C), the Min system appears to remain in definite eigenmodes. That is, the spatial patterning of MinD (y-direction in the kymogram) can always be approximated by a simple cosine wave, ${f_n}\left( x \right) = \cos \left( {\frac{{n\pi }}{l}x} \right)$ where n is the spatial mode number and l is the cell length.

Figure 4A shows a plot of the dominant spatial mode of the system calculated using full simulations of the non-linear equations as the length of the cell is increased (shown in red). The cell length at which transitions occur lags that determined via linear stability analysis (shown in blue). However, the patterning of the system still transitions through the harmonic modes of the cell geometry as the length of the cell increases. That is, the first order mode transitions to the second which transitions to the third and so on.

If the dispersion curve around the maximum is approximately symmetric, then the distance between critical transition lengths in one dimension should be approximately constant. The local symmetry is demonstrated in Figure 2, where the whole dispersion curve is highly asymmetric, however the area near the maximum, that is the region between the first two non-zero blue points is approximately symmetric. When the cell length is such that two consecutive modes are equally destabilised (equal values of ${\lambda _{\max }}\left( k \right)$), that is, at a mode transition, the two corresponding values, k_n and k_n+1 are equally spaced about the maximum of the dispersion curve. This separation in the wavenumber dimension, k, is maximal for the transition between the first two modes (between $k = 0.94\mu {m^{ - 1}}$ and $1.88\mu {m^{ - 1}}$ in Figure 2). If the dispersion relation is symmetric in this region, then the cell lengths at which mode transitions occur will be equally spaced.

Full simulations of the non-linear equations describing the Min model show that the transitions are approximately equally spaced (Figure 4A). This makes it possible to approximately describe the mode transitions of the Min system by fitting for the maximum value of a heuristic dispersion relation. Linearly fitting the mode number plot from Figure 4A, gives a maximum amplification at a wavenumber, k, of $0.78\mu {m^{ - 1}}$. From this value the predicted mode transitions are shown in Figure 4B. An overlay of the mode transitions obtained from the full simulations with those predicted from the heuristic dispersion relation are shown in Figure 4C. Yellow regions indicate where both are identical.

Thus, it is possible to obtain a heuristic dispersion relation from models or experimental data on the mode transitions of a system where the cell geometry resembles a one-dimensional object (for example, a capsule shaped cell, such as wild type E. coli). Based on our observations obtained from the full simulations of the non-linear equations, this heuristic dispersion relation will have the following properties: (i) ${\lambda _{\max }}\left( k \right)$ becomes positive for $k < 3.35\mu {m^{ - 1}}$, thus, no patterning is seen for cell lengths less than $0.935\mu m$; (ii) for k values between $1.16\mu {m^{ - 1}}$ and $3.35\mu {m^{ - 1}}$ (a cell length between $0.935\mu m$ and $2.7\mu m$), ${\lambda _{\max }}\left( k \right)$ is real, thus, a stationary MinD pattern will be observed that corresponds to the first order mode; (iii) for wavenumbers less than $1.16\mu {m^{ - 1}}$ (a length greater than $2.7\mu m$) ${\lambda _{\max }}\left( k \right)$ is complex, thus, the system will oscillate in the mode whose wavenumber closest to the maximum in the dispersion relation $k = 0.78\mu {m^{ - 1}}$. This heuristic dispersion relation is shown in Figure 5.

3.2.6. Testing the heuristic dispersion relation on an arbitrary shape

The coordinate free approach to linear stability analysis shows that the dispersion relation is independent of geometry. Thus, it should be possible to utilise a heuristic dispersion relation, such as the one derived above, for cells of an arbitrary shape.

To test this, full simulations of the non-linear equations describing the Min system were run on an arbitrary domain (two-dimensional cell shape) and compared to the predictions based on the heuristic dispersion relation. The arbitrary shape was created by randomly generating a Fourier basis expansion for a parametric outline (see Methods for details). The shape that resulted from this process can be seen in Figure 6. In each panel, maximum MinD concentration is shown as yellow, with minimum MinD as blue. For oscillating patterns, green regions represent nodes where the MinD concentration does not vary.

The size of the domain is controlled by linearly scaling the shape in all directions using the scaling factor C. It is worth noting that allowed wavenumbers k (related to the eigenvalues, -k², of the Laplace-Beltrami operator), scale inversely proportional to the scaling factor (C). Halving the size of the shape, that is halving C, will double all wavenumbers and vice versa. The lowest wavenumber for the shape is $k = 1.71\mu {m^{ - 1}}$ when no scaling is performed (C = 1).

According to the predictions based on the heuristic dispersion relation, scaling the shape to 51% (C = 0.51) of its original size should increase the lowest non-zero wavenumber describing the shape to the critical value of $k = 3.35\mu {m^{ - 1}}$. This corresponds to the minimum size shape in which patterning of MinD occurs. As shown in Figure 6A, by simulating the full non-linear set of equations describing the Min system shapes that are smaller than C = 0.45 display no patterning, while shapes that are larger than C = 0.55 pattern in a stationary first order mode.

Similarly, a scale factor of C = 1.47 will decrease the lowest wavenumber of the shape to the critical value of $k = 1.16\mu {m^{ - 1}}$, below which the stationary pattern begins to oscillate. In Figure 6B full simulations of the non-linear equations show that cells smaller than this maintain a stationary pattern where maxima and minima remain permanently in a single location. For example in Figure 6Bi with a scale factor of C = 1.4, high concentration always remains on the right side of the cell. The concentration of MinD in cells larger than this oscillates in time; see Figure 6Bii where a scale factor of C = 1.5 causes the high concentration to oscillate back and forth.

Figure 6C shows the MinD distribution as the shape grows via scaling C from 2.2 to 7.98. For each pair (upper and lower), the left and right panels represent the extremes of the MinD distribution as over one period of oscillation. The upper panels, labelled “Prediction”, are the first six eigenmodes, that is, eigenfunctions of the Laplace-Beltrami operator, obtained from the heuristic dispersion relation while the lower panels are the results of full simulation of the non-linear equations. For C = 2.2 (Figure 6C), the heuristic dispersion relation selects the first eigenmode (upper panels) which accurately represents the result of the full simulation (lower panels). For successive eigenmodes (C = 3.28, 4.48 etc), correspondence between the MinD distribution from the full simulation and the eigenfunction selected by the heuristic dispersion relation is not as accurate as for the first eigenmode (Figure 6C). However, decomposition of the MinD distribution obtained from the full simulation into eigenfunctions shows that the eigenfunction selected by the heuristic dispersion relation is the dominant term (Figure 7). This is seen as the maxima, in yellow, and minima, in blue, seen in the prediction (top row), are always present in the respective simulation (bottom row). However, areas where harmonics have large nodal regions (regions of green) are not seen in simulations (for example in the second mode of Figure 6C). Instead these areas pattern to form unexpected maxima and minima in the respective simulations. These additional maxima and minima are coming from contributions due to neighbouring harmonics.

The transition through modes as the size of the shape is increase is shown in Figure 7. For the full simulations, the MinD patterns were decomposed into eigenfunctions to determine the dominant mode. The dominant mode of patterning for full simulations at each size is shown in Figure 7A (note that full simulations were run for scale factors incremented by steps of $\Delta$C = 0.25, hence, transitions are seen as lines with a finite slope). In comparison the heuristic prediction of mode transitions assuming mode closest to the wavenumber of $k = 0.78\mu {m^{ - 1}}$ is dominant is shown in Figure 7B. The overlay of the transitions calculated using transitions (red), with the heuristic prediction (green), and areas where both are identical (yellow) is shown in Figure 7C. The smaller step size for predictions versus simulations leads to the lines between points being steeper in the prediction plot than for the simulations. With this in mind, both curves follow a similar trend, suggesting the strong predictive power of this heuristic approach.

3.3. Application to experimental results

The computer simulations above show that a full simulation of the non-linear equations describing the Min system can be used to generate a heuristic dispersion relation. The simulations indicate that a separation between the reaction terms and the geometry of the domain (cell) is maintained. Furthermore, the pattern observed in the full simulations will be dominated by the harmonic mode (eigenfunction of the Laplace-Beltrami operator) whose eigenvalue, -k², corresponds to the wavenumber, k, closest to the positive peak in the heuristic dispersion relation. This approximation is particularly good for the low order modes. It is important to determine whether this approach can be applied to real experimental data to predict MinD patterning. In particular, we wish to test this against data where the shape of the cell has been modified.

3.3.1. Rectangle harmonics approximate Min oscillations in rectangular cells

Recently, Wu et al. studied the patterning of the Min system in E. coli cells which they were able to deform into rectangles using PDMS micro-chambers [23]. As a part of this study the length and width of the rectangles were varied to investigate the effect of geometry on the patterning of the Min system. Wu et al. identified seven unique types of patterning for E. coli with rectangular shapes [23]. Examples of each of these patterning types are shown in the theory column of Figure 8.

The seven identified types of patterning can all be closely approximated by either individual low order harmonics of a rectangle or the sum of two sequential harmonics (single modes and mixed modes, respectively, Figure 8). The harmonic modes of a rectangle can be determined analytically and are of the form ${f_{n, m}}\left( {x, y} \right) = \cos \left( {\frac{{n\pi }}{l}x} \right)\cos \left( {\frac{{m\pi }}{w}x} \right)$ where each mode is indexed by two integers (n, m) corresponding to the x- and y-directions. These mode numbers are shown in the left hand column of the Figure 8. The pure modes ((1, 0), (0, 1), (2, 0) and (1, 1)) are shown in the upper part of Figure 8, while the mixed modes are shown in the lower part of the Figure.

For the pure modes, the experimental observations and the theoretical predictions show a close correspondence (Figure 8). In the experiments [23], patterning that we associate with mixed modes occurs in degenerate shapes. That is, shapes where the x- and y-dimensions are such that two different modes are simultaneously amplified in a way that neither dominates the other (examples include identical x- and y-dimension and rational multiples of one dimension). The first pattern shown in the mixed mode section of Figure 8 is generated by the summation of the (1, 0) mode (shown on the first line of the single mode section) and a (0, 1) mode (shown on the second line of the single mode section). These two modes are added in phase. That is, the modes are synchronized so that when the (1, 0) mode has a maximum at the left hand edge of the rectangle, the (0, 1) mode has it maximum value at the bottom edge. If these two modes are added out of phase, that is, when the (1, 0) has its maximum at one end, the (0, 1) mode is half way between ends, the resulting patterning has a maxima that rotates around the edge of the rectangle. This rotational patterning is shown on the second line of the mixed mode section in Figure 8. Both of these modes exist when the rectangle is essentially a square so that the two individual modes have the same wavenumber, k.

In rectangles which are twice as long as wide, the second order mode in the x direction ((2, 0) shown on third line of the individual modes) will have the same wavenumber, k, as a first order mode in the y direction ((0, 1) shown on the second line of individual modes). The result of their mode mixing is shown on the third line of the mixed mode section in Figure 8.

3.3.2. Patterning of the Min system in deformed E. coli is determined by cell shape

The work of Männik et al. squeezed E. coli cells into narrow nano-fabricated channels where they exhibit highly irregular shapes and large volumes [24]. As a part of their studies of the robustness of cell division machinery in these deformed cells, they collected time-lapse data of the Min protein system.

To compare the time-lapse measurements to the heuristic dispersion relation analysis, the outlines of the cells reported in Männik et al. [24] were traced and their respective harmonics calculated. To account for time averaging across the entire oscillation, the absolute value of the relevant harmonic was taken so that both maxima and minima appear red at the same time. The results are shown in Figure 9.

The left hand cell in Figure 9 displays a simple oscillation (Figure 9A), so the prediction of the oscillation fits the time-lapse measurement very well (Figure 9C). In the right hand image (Figure 9B), the Min system is in a higher, more complicated mode. As a result, the prediction is less accurate but still confers the main feature of the patterning (Figure 9D).

4. Discussion

In all domains of life, it is clear that the accurate location of subcellular structures, including individual macromolecules, is essential for cellular function. In most instances, the mechanisms that determine subcellular localisation are poorly characterised. A prime example of intracellular localisation is the cell division machinery. The complex systems that underpin cell division couple multiple protein-protein and protein-membrane interactions to the geometry of a cell in order to precisely determine the position and timing of cell division. This precision is required for successful reproduction of organisms.

A full characterisation of the systems responsible for positioning the cell division machinery is lacking in most instances, with the possible exception of the E. coli Min system, where several mathematical models have been proposed [16,17,18,19,20]. The Min protein system establishes an oscillating spatial pattern that is important in positioning the FtsZ division ring. The Min protein system can be described by reaction-diffusion equations that produce spatiotemporal patterns [1,2]. Thus, it is important to develop heuristic tools that can provide a means of predicting the behaviour of the patterning systems similar to Min based on models or experimental data.

Taking a coordinate-free approach to linear stability analysis provides a method to separate reaction terms from cell geometry. This separation implies that it is possible to predict how a system will pattern in any shape by either deriving the dispersion relation using linear stability analysis, or heuristically, by knowing how the system patterns in a single family of cell shapes. Predictions are in the form of harmonic modes determined by the geometry of the cell shape.

Using a full mathematical model of the Min system [20], linear stability correctly predicts the onset of patterning but fails to predict temporal oscillations and the cell lengths at which transitions between successive modes occur. Although linear stability analysis fails in a predictive sense, the underlying tenets of linear stability analysis still hold. In particular these are: that the patterning of MinD approximates one of the harmonic modes of the cell shape; that there is a sequential transition through harmonic modes with increasing cell size; and finally the underlying reactions can be represented by a heuristic dispersion relation, separating the reactions from the geometry of the cell.

Using the results of a full simulation of the non-linear equations that describe the Min system on a one-dimensional model cell, we have derived a heuristic dispersion relation. Using this heuristic dispersion relationship, we are able to closely approximate the results of the full simulation of the non-linear equations. Furthermore, we have been able to use this heuristic dispersion relation to predict the Min protein patterns that have been experimentally observed for E. coli cells that have been deformed into either rectangles or arbitrary planar shapes.

Based on this, it should be possible to construct heuristic dispersion relations for any cell patterning system based on either a set of mathematical simulations or experimental data showing the behaviour of the patterning proteins during cell growth. Utilising the separation of reaction terms from cell geometry, it should then be possible to use this heuristic dispersion relation to predict patterning on an arbitrary cell shape.

The discovery of new examples where protein patterning is consistent with a heuristic dispersion relation has implications for the properties of the underlying protein system. Such a system is likely to involve a Turing instability and must be actively driven in order to dynamically maintain patterning. Thus, the underlying the protein system must contain an energy source, most probably an ATPase or GTPase. Given the requirement for non-linear reaction dynamics to produce patterning, it is also probable that the system will contain proteins capable of some degree of oligomerisation.

It is important to appreciate the limitations of applying reaction-diffusion formalism to intracellular localisation problems. The computations in this work have been performed on cells of fixed size that are assumed to start with random initial protein distributions. In reality, cells are growing and have existing patterns that are likely to be changing at different rates. This has the potential to introduce hysteresis effects into the determination of patterning seen in a system, most notably in the form of mode doubling [34]. In real cells, there is likely to be interplay between hysteresis, which will cause the biasing of theoretically unfavourable modes, and stochastic effects, which will degrade hysteresis over time and push the system back to the most favourable pattern.

5. Conclusion

The core tenets of linear stability analysis utilized in this paper require no knowledge of the underlying protein dynamics, and so can potentially be used as an assay for the existence of unidentified Turing patterning systems. While there is no guarantee other intracellular Turing patterns would exhibit these characteristics (the existence of a heuristic dispersion relation that accurately predicts patterning), given the fundamental nature of the derivation, it seems probable that this will be the case.

Acknowledgments

JCW is a recipient of an Australian Postgraduate Award.

Conflict of Interest

The authors declare that they have no conflicts of interest with the contents of this article.