## Abstract

The distributions of many proteins in rod-shaped bacteria are far from homogeneous. Often they accumulate at the cell poles or in the cell centre. At the same time, the copy number of proteins in a single cell is relatively small making the patterns noisy. To explore limits to protein patterns due to molecular noise, we studied a generic mechanism for spontaneous polar protein assemblies in rod-shaped bacteria, which are based on cooperative binding of proteins to the cytoplasmic membrane. For mono-polar assemblies, we find that the switching time between the two poles increases exponentially with the cell length and with the protein number. This feature could be beneficial to organelle maintenance in ageing bacteria.

## 1. Introduction

Although usually small in size, bacteria display a remarkable degree of spatial order [1,2]. In the simplest case, proteins are in a controlled way unevenly distributed between the two halves of rod-shaped bacteria like *Escherichia coli* or *Bacillus subtilis*. Examples for asymmetric distributions are provided by the accumulation of chemotactic receptors at the poles of *E. coli* [3], the localization of Spo0J/Soj proteins close to the poles of *B. subtilis* [4,5], or the Min proteins in *E. coli*, which periodically switch between the two cell halves [6]. More refined structures are the Z-ring that localizes at the site of cell division [7], the localization as a string of pearls of magnetosomes involved in bacterial magnetotaxis [8] and the distribution of flagella of various motile bacteria [9].

Typically, the spatial distribution of proteins is linked to their function. Tom Duke was among the first to realize the importance of spatial structures for signalling. In the context of bacterial chemotaxis, he argued that the formation of receptor complexes increases sensitivity [10]. Another example is provided by the Min protein oscillations that are used to select the site of cell division [11]. To give just one more example, the linear arrangement of magnetosomes produces a strong enough magnetic dipole moment for bacteria to orient along the field lines of the geomagnetic field [8].

As a consequence of their functional importance, these spatial structures need to be maintained for a proper working of the cell. However, usually, the relatively small copy number of proteins in a cell, which typically ranges between less than 10 to a few thousands, implies that molecular noise may severely restrict their lifetime. For some structures, noise does not seem to interfere with their function. For example, the Spo0J/Soj proteins stochastically switch between the two nucleoids without impairing the fitness of *B. subtilis* [4,5]. In other cases, molecular mechanisms ensure the stability of protein patterns in spite of small copy numbers. For instance, numerical analysis of the Min protein dynamics has revealed that the oscillatory pattern is robust against molecular noise [12–14].

The influence of molecular noise on bacterial processes has intensively been studied during recent years [15]. In this context, a particular focus has been put on bistable systems [16], although spatially extended systems have received less attention. One exception is provided by the Spo0J/Soj system mentioned above [4,5,17]. Also the Min proteins in *E. coli* were studied in this context and found to stochastically switch between the two cell halves in sufficiently small cells after moderate over-expression [18–20]. So far, however, we lack a general understanding of how noise affects spatial protein distributions and theoretical tools for their analysis still need to be developed further.

The spatial structures mentioned above have in common that the proteins in question assemble on a support, for example, the membrane or the nucleoid. Consequently, spatial cues on these scaffolds might underlie the formation of the protein aggregates. In an extreme case, the proteins would not interact with each other but rather move in a potential landscape imposed by spatial cues on the scaffold. In such cases, the stability of the pattern can be analysed with standard Kramers theory [21]. Furthermore, entropic effects can lead to polar localization, for example, through nucleoid occlusion [22]. An alternative concept explaining structure formation is protein self-organization [2]. Mechanisms for self-organization have been proposed, for example, for Spo0J/Soj [17] and the Min proteins [23]. For the latter, the possibility of self-organization has been demonstrated in reconstitution experiments [24–26]. Currently, there is no general framework for analysing the impact of noise on spatial self-organization.

As a particular example for the impact of noise on cellular protein patterns, we study in this work the switching between two self-organized states of a spatially extended system in the weak-noise limit. In this case, typically, one cannot use Kramers rate theory, which relies on the existence of a potential [21]. Instead, a generalization of Kramers' theory based on a pseudo-potential method can be employed [27–30]. In a biological context, this method has been applied to bistable genetic switches [31,32] and to bidirectional transport of molecular motors [33,34]. Our study is motivated by the heterogeneous protein distributions in rod-shaped bacteria described above. Heterogeneity of the distributions results from cooperative protein attachment to the membrane, which is motivated by studies of the Spo0J/Soj dynamics [17] and the Min system [35]. The molecular origin of cooperative binding remains to be understood, but assuming this process has been very successful in these contexts. After defining the model, we will first perform stochastic simulations. To get further insight, we will perform a mean-field analysis and then follow the approach by Hildebrand & Mikhailov [36] to establish the corresponding Fokker–Planck equation for the dominant modes. Employing Wentzell–Freidlin theory, we will solve the Fokker–Planck equation and obtain the switching time as in [29,30]. The work concludes with some remarks about possible functional implications and generalizations.

## 2. Stochastic dynamics for molecules forming membrane clusters

### 2.1. The chemical master equation

In the following, we consider the following processes: (i) cytoplasmic molecules can bind to the membrane, with binding being favoured in regions where membrane-bound molecules are already present, (ii) spontaneous detachment of membrane-bound molecules, and (iii) diffusion of cytoplasmic and membrane-bound molecules (figure 1). In the systems mentioned in the Introduction, bound molecules are released by the action of an antagonist. Including the action of an antagonist, however, leads to quite involved dynamics. We will thus focus here on the simpler case of spontaneous unbinding and leave the case of induced unbinding to a subsequent publication. Furthermore, we will reduce the dynamics to occur along a line, which we expect to capture the essence of the protein dynamics along the long axis of a rod-shaped bacterium as long as the protein patterns are invariant with respect to rotations around the axis.

The system is of length *L* and contains *Ω* particles. The number of particles in our system can be compared with the number of molecules *N*_{tot} in a bacterium. Assuming that, in the bacterium, the distribution is homogeneous with respect to rotations around the bacterium's long axis and assuming the same total density in the azimuthal as in the longitudinal direction, we get , where *d* is the bacterial diameter. For example, for typical size of *E. coli* with and , , such that for a concentration of the order of a few micromolar, we have *Ω* of the order of 50–100. In this work, we will assume *L* and *Ω* to be constant during one simulation. The effects of bacterial growth and the production of proteins will be studied in later work.

For our calculations, we decompose the system into compartments of length *ℓ* and assume that particles within a compartment are well mixed. The rate of spontaneous attachment to the membrane is denoted by *ω*_{a}, whereas the rate of spontaneous detachment is *ω*_{d}. The rate of cooperative membrane-binding in compartment *i* is given by *ω*_{c}*n*_{m,i}(*n*_{m,i} − 1)/(*Ω*^{2}*ℓ*^{2}). In this expression, *n*_{m,i} is the number of membrane-bound molecules in compartment *i*; the corresponding number of cytosolic molecules will be denoted by *n*_{c,i}. We make the scaling of the cooperative binding term with the compartment size and the total molecule number explicit, because we will later consider the limit of large molecule numbers *Ω* → ∞ and the continuum limit *ℓ* → 0. Note that in order to compare systems with the same density at different lengths in the limit *Ω* → ∞, we will scale *ω*_{c} as *L*^{2} in §§3 and 4.

The corresponding chemical master equation for the evolution of the joint probability *P*({*n*_{c,i}}_{i}_{=1,…,N}, {*n*_{m,i}}_{i}_{=1,…,N}; *t*) can be written as
2.1In this expression, denote the particle creation and annihilation operators in compartment *i* such that, for example, . The corresponding operators for membrane-bound compartments are denoted by . Products of and indicate that a particle is taken away from one compartment or state and put into another compartment or state, respectively. The terms in the first three lines of the chemical master equation account for attachment to and detachment from the membrane, the next lines capture diffusion of cytoplasmic and membrane-bound molecules. The attachment rates for compartment *i* are proportional to the number of cytosolic molecules *n*_{c,i}, the corresponding detachment rate proportional to the number of membrane-bound molecules *n*_{m,i}, and similarly for the ‘ hopping rates’ in the diffusion terms.

### 2.2. Numerical solution of the chemical master equation

We numerically analyse the chemical master equation by combining the Gillespie algorithm with the next subvolume method [37]. Explicitly, in a given state, we determine for each compartment *j* the total rate of a change in the number of either the bound or unbound particles. For each compartment, we then draw a random number to determine the time at which the next event occurs and perform the corresponding action for the compartment with the shortest waiting time. We then go back and update the total rates of all compartments and continue as before.

In figure 2, we present the results of a simulation for a system length of 5.2 µm. The kymograph in figure 2*a* shows that the particles accumulate in the vicinity of one end or ‘cell pole’. At random times, the particles switch to the opposite pole. The distribution of switching times has a mean value and is rather broad (figure 2*c*) with a standard deviation *σ*_{τ} = 5487 s. The time-averaged profile of particles in one cell half agrees well with the profile obtained from the mean-field equations (3.14) and (3.15) (figure 2*b*) which will be introduced and discussed later. It should be noted that, as for all continuum theories, the agreement between the mean-field theory and the stochastic simulations depends on the choice of *ℓ* in the stochastic simulations. Each bin should contain on average a sufficiently large number of molecules such that a continuum approach is appropriate. At the same time, *ℓ* should be small compared with spatial features of the protein distribution.

## 3. The functional Fokker–Planck equation

In order to go beyond the numerical analysis of the system, we will derive a functional Fokker–Planck equation from the chemical master equation (2.1) following essentially the method of Hildebrand & Mikhailov [36]. Explicitly, we will perform in a first step an expansion in terms of the inverse molecule number 1/*Ω*, which is a form of van Kampen's system size expansion. In a second step, we will take the continuum limit *ℓ* → 0. The resulting equation will first be analysed in the deterministic limit *Ω* → ∞. Then, we will use a Galerkin ansatz to obtain an equation that is amenable to analysis of the influence of molecular noise on the system behaviour, which will be carried out in the next section.

### 3.1. System size expansion and continuum limit

Let us consider the case of . In each compartment with large occupancies, and , we can approximate the creation and annihilation operators of cytosolic and membrane-bound molecules, and , respectively, by 3.1and 3.2

Formally, this corresponds to an expansion in the inverse molecule number up to second order. For large but finite values of *Ω*, the above condition might not be fulfilled for all compartments at all times and higher orders of the approximation should in principle be taken into account. Still, we will see that the ensuing equations give a quantitative account of the system behaviour. In the same spirit, we will also replace from now on the factor *n*_{m,i} − 1 in the cooperative binding term by *n*_{m,i}. Explicitly, the Fokker–Planck equation resulting from the chemical master equation (2.1) reads as follows:
3.3

We will not use the Fokker–Planck Equation in this form, but rather go on by taking the continuum limit *ℓ* → 0. To this end, we introduce the densities *c _{i}* =

*n*

_{c,i}/(

*Ωℓ*) and

*m*=

_{i}*n*

_{m,i}/(

*Ωℓ*) for

*i*= 1, … ,

*N*. They give the fraction of all molecules that are cytosolic or membrane-bound in compartment

*i*, such that . Note that

*c*like

_{i}*m*is a line-density, because we assume uniformity of the distribution in radial direction of the cell.

_{i}In the limit *ℓ* → 0, the densities can be replaced by continuous functions *c* and *m* with *c*(*iℓ*) = *c _{i}* and

*m*(

*iℓ*) =

*m*. Correspondingly, sums in equation (3.3) are replaced by integrals, , and partial derivatives by functional derivatives according to 3.4and 3.5and analogously for the other partial derivatives.

_{i}The resulting functional Fokker–Planck equation for the probability functional *P*[*c*(*x*), *m*(*x*), *t*] reads as follows:
with
3.7
3.8
3.9
3.10
3.11where
3.12and
3.13The terms containing *A*_{c} and *A*_{m} describe the ‘convective’ or deterministic part, whereas the ‘diffusion’ terms containing *B*_{cc}, *B*_{cm}, *B*_{mc} and *B*_{mm} account for the influence of noise. Note that the latter are proportional to the inverse total molecule number 1/*Ω*. In the limit *Ω* → ∞, we get *P*[*c*(*x*), *m*(*x*), *t*] = *P*_{0}[*c*(*x*, *t*), *m*(*x*, *t*)], where *c*(*x*, *t*) and *m*(*x*, *t*) obey the time evolution equations in the deterministic limit *∂ _{t}c* =

*A*

_{c}and

*∂*=

_{t}m*A*

_{m}, and where

*P*

_{0}[

*c*(

*x*),

*m*(

*x*)] =

*P*[

*c*(

*x*),

*m*(

*x*),

*t*= 0].

### 3.2. The deterministic limit

We start the analysis of the functional Fokker–Planck equation by considering the deterministic limit. Explicitly, the dynamic equations for the particle densities *c* and *m* read as follows:
3.14and
3.15They are complemented by no-flux boundary conditions on the diffusion currents −*D*_{c}*∂ _{x}c* = 0 and −

*D*

_{m}

*∂*= 0 at

_{x}m*x*= 0 and

*x*=

*L*.

Introducing the dimensionless time *t’* = *tω*_{a} and space *x’* = *x*/*λ* with *λ*^{2} = *D*_{c}/*ω*_{a} and dimensionless densities *c’* and *m’* through and , the dynamic equations can be written in the dimensionless form
3.16and
3.17where we have dropped the primes for ease of notation. The dimensionless diffusion constant is *D* = *D*_{m}/*D*_{c} and the dimensionless detachment rate *k* = *ω*_{d}/*ω*_{a}. Note that also the system size *L* is now measured in units of *λ*.

As mentioned in §2, we choose from now on to compare systems with the same particle density. The deterministic equations have a stationary homogeneous solution *c* = *C*_{0} = const. and *m* = *M*_{0} = const. that are then determined by and
3.18For the parameter values used in this work, there is only one real solution to this equation.

For a linear stability analysis of the homogeneous stationary solution, we expand the densities
3.19and
3.20and keep only linear terms in *C _{n}* and

*M*,

_{n}*n*= 1, 2, … in the dynamic equations 3.21where

*q*=

_{n}*nπ*/

*L*and

*n*= 1, 2, …

The homogeneous state is stable for *k* ≤ *k*_{c,1} or *k* ≥ *k*_{c,2} with some critical values *k*_{c,1} and *k*_{c,2}. In between these two values, there is a critical system length *L*_{c}(*k*) beyond which the homogeneous distribution becomes unstable through a pitchfork bifurcation (figure 3*a*). The first mode to become unstable is the one with *n* = 1 corresponding to situations in which the molecules accumulate at a pole. As the system is invariant under the transformation *x* → *L* − *x*, there coexist two mirror-symmetric solutions. One of them is chosen by spontaneous symmetry breaking (figure 3*b*).

Increasing the system length beyond *L*_{c}, there is a second bifurcation at *L* = *L*_{c,2}(*k*) through which two new solutions appear. They correspond to states where the proteins either pile up in the cell centre or accumulate at both poles. Numerical integration of the deterministic equations (3.16) and (3.17) indicates that only the state with proteins clustering in the centre is stable (figure 3*b*).

We can make a Galerkin ansatz and truncate the series (3.19) and (3.20) at a finite value of *n*. The dynamic equations resulting after insertion of the truncated series in equations (3.16) and (3.17) are given in appendix A. It turns out that for a series with *n* = 2 gives a faithful description of the system dynamics (figure 3*c*). Furthermore, the modes *C*_{0}, *C*_{1}, *C*_{2} and *M*_{0} relax on faster time scales than the modes *M*_{1} and *M*_{2}, such that we can make an adiabatic approximation and express the former in terms of the latter. The dynamic equations thus only depend on *M*_{1} and *M*_{2}.

In figure 4, we present the flow field for the dynamic equations for *M*_{1} and *M*_{2} after adiabatic elimination of *C*_{0}, *C*_{1}, *C*_{2} and *M*_{0}. For *L*_{c} ≤ *L* ≤ *L*_{c,2}, the system has two symmetric stable and one hyperbolic fixed point (figure 4*a*). The hyperbolic fixed point is at (0,0,0,0) and corresponds to the homogeneous state. The symmetric fixed points have such that they correspond to protein accumulations at either of the two poles. As soon as *L* > *L*_{c,2} two new fixed points appear. These hyperbolic fixed points have *M*_{1} = 0 and correspond to symmetric protein distributions with, respectively, accumulation in the centre and simultaneous accumulation at both poles. Further increase in *L* leads to two new hyperbolic fixed points. At the same time, the fixed point with *M*_{1} < 0 becomes stable. It should be noted, though, that the corresponding symmetric distribution is not observed in the stochastic simulations, indicating that this fixed point remains hyperbolic if higher modes are taken into account. The two new hyperbolic fixed points are connected to the stable fixed points by heteroclinic orbits (figure 4*b*).

### 3.3. Galerkin approximation for the Fokker–Planck equation

We will now exploit the findings of the previous subsection and use the Galerkin approximation for the Fokker–Planck equation (3.6). We will use the same dimensionless parameters and densities as above. Employing equations (3.19) and (3.20), we can transform the functional *P*[*c*(*x*), *m*(*x*); *t*] into a function of the coefficients *C _{n}* and

*M*,

_{n}*n*= 0, 1, 2, … . The corresponding Fokker–Planck equation can be derived from equation (3.6). To this end, we first express the functional derivatives with respect to

*c*(

*x*) and

*m*(

*x*) in terms of partial derivatives with respect to

*C*and

_{n}*M*,

_{n}*n*= 0, 1, 2, … . Considering a variation

*δ C*of mode

_{n}*n*of cytosolic molecules yields

*δc*(

*x*) =

*δ*C

_{n}cos(

*nπx*/

*L*). If denotes the function of

*C*and

_{n}*M*corresponding to a functional

_{n}*F*[

*c*(

*x*),

*m*(

*x*)], then 3.22and analogously for the derivatives . Note that variations of

*C*

_{0}and

*M*

_{0}are not independent as

*C*

_{0}+

*M*

_{0}= const. From the above expression it follows that 3.23and 3.24and similarly 3.25 3.26 3.27

Inserting expansions (3.19) and (3.20) into equation (3.6), using the expressions for the functional derivatives just derived, and performing the integrals, we obtain a Fokker–Planck equation for the probability density *P*({*C _{n}*}

_{n}_{=0,1,…},{

*M*}

_{n}

_{n}_{=0,1,…};

*t*), where again

*C*

_{0}+

*M*

_{0}= const. It is of the form 3.28We will now apply the Galerkin approximation and truncate the series at

*n*,

*m*= 2. Furthermore, we will make the adiabatic approximation and express

*C*

_{0},

*C*

_{1},

*C*

_{2}and

*M*

_{0}in terms of

*M*

_{1}and

*M*

_{2}by using the deterministic equations (A 3)–(A 6) in steady state. Eventually, we arrive at 3.29where the probability distribution

*P*≡

*P*(

*M*

_{1},

*M*

_{2};

*t*) now depends only on

*M*

_{1}and

*M*

_{2}and where

*∂*≡

_{k}*∂*/∂

*M*with

_{k}*k*= 1, 2. Furthermore, and summation over repeating indices is understood. The expressions for the drift velocity are those given in the deterministic limit (equations (A 7) and (A 8)). The expressions for the diffusion matrix are given in appendix A. Note that the elements

*D*depend on

_{kl}*M*

_{1}and

*M*

_{2}.

In the next section, we will use the Fokker–Planck equation (3.29) to determine the average switching time.

## 4. Estimation of the switching time in the weak-noise limit

The switching time can equivalently be interpreted as the mean first-passage time (MFPT) of a multi-dimensional escape problem. However, in contrast to most cases studied in the literature, the deterministic dynamics is in the present case not determined by a potential, so that we cannot use Kramers (or Eyring) rate theory; instead, we will employ Wentzell–Freidlin theory [29,30]. In a first step, we will present the necessary expressions and then compare their solution to data obtained from stochastic simulations. We continue to use the dimensionless parameters introduced above.

### 4.1. Probability flux across the separatrix

As we have seen above, in the deterministic case there is a region in parameter space where two stable mirror-symmetric distributions coexist. The two basins of attraction are separated by a separatrix along *M*_{1} = 0, which contains another fixed point, namely the homogeneous state *M*_{1} = *M*_{2} = 0. It is a hyperbolic fixed point for which the stable manifold coincides with the separatrix. The aim is now to calculate the MFPT *τ* in the limit of weak noise.

To this end, we consider the Fokker–Planck equation (3.29) with absorbing boundary conditions along the separatrix. We will calculate the total probability current across the separatrix associated with the eigenfunction *P*_{1} of the Fokker–Planck operator that decays most slowly. Normalizing this current to the total probability to find a particle on one side of the separatrix yields
4.1where summation with respect to the index *k* is understood. Note that *τ* is just the characteristic time on which the eigenfunction *P*_{1} relaxes, that is, it is the inverse of the corresponding eigenvalue. In lowest order in 1/*Ω*, the eigenfunction *P*_{1} can be replaced by the solution *P*_{s} to the steady-state Fokker–Planck equation [38].

The steady state of equation (3.29) cannot be calculated exactly. Instead, we will make a Wentzel–Kramers–Brillouin (WKB) ansatz and write the steady-state probability distribution as
4.2with the classical action (or quasi-potential) *S*. The equation for the action is obtained by inserting the WKB ansatz into equation (3.29) and by considering the terms of first order in 1/*Ω*. We get an equation of the form of a Hamilton–Jacobi equation *H*(*M*_{1}, *M*_{2}, ∂_{1}*S*, *∂*_{2}*S*) = 0 with the Wentzell–Freidlin Hamiltonian
4.3where *p _{k}* =

*∂*is the momentum conjugated to

_{k}S*M*,

_{k}*k*= 1, 2. The action is then obtained from solving first the canonical dynamic equations, which in our case read 4.4and 4.5and then the dynamic equation 4.6It remains to determine the prefactor

*K*in equation (4.2), which can be obtained from the terms of zeroth order in 1/

*Ω*in equation (3.29). Explicitly, 4.7

In this expression, denotes the Hessian of the action *S*, that is, *W _{kl}* =

*∂*. It can be determined without knowledge of the action by solving its time evolution equation, which follows a (generalized) Riccati equation 4.8

_{k}∂_{l}SIn figure 5, we present steady-state probability distributions for various system lengths. The probability is maximal in the vicinity of the stable fixed points. In the hyperbolic fixed point at (0,0,0,0), the probability distribution has a saddle point (slight deviations are due to interpolation errors) and the heteroclinic orbit indicated in yellow follows the ridge of probability distribution.

We can exploit the form of the probability distribution to further simplify equation (4.1). The integrand in the denominator of the above expression for *τ* is obviously dominated by the region around the stable fixed point, whereas the integrand in the numerator is dominated by the region around the hyperbolic fixed point *M*_{1} = *M*_{2} = *p*_{1} = *p*_{2} = 0. In these regions, we can make a Gaussian approximation of the probability density or, equivalently, use a quadratic approximation of the Hamiltonian (4.3):
4.9where *λ _{kl}* ≡

*∂*and the coefficients are evaluated in the hyperbolic or stable fixed points of equations (A 7) and (A 8), respectively. With this approximation, the switching time is estimated as 4.10In this expression, the quantities with a superscript ‘hyp’ are evaluated in the hyperbolic fixed point, whereas a superscript ‘st’ indicates evaluation in the stable fixed point. Note that if

_{k}u_{l}*ω*

_{c}scales with

*L*

^{2}, then the prefactor of

*S*

^{hyp}in the exponential is independent of

*L*.

Finally, the expression for and can be obtained from equation (4.8) by setting the time-derivative equal to zero and using *p*_{1} = *p*_{2} = 0. It yields
4.11The elements *D _{mn}* and the partial derivatives

*∂*again have to be evaluated at the respective fixed point.

_{m}u_{n}We have now given all the elements necessary for estimating the switching time *τ*.

### 4.2. Comparison of the estimated switching time with stochastic simulations

According to equation (4.1), we need to evaluate various quantities in the fixed points of equations (A 7) and (A 8). In the hyperbolic fixed point, the approximate Hamiltonian (4.9) can be obtained analytically. The matrices *λ* and are diagonal with non-zero elements
4.12
4.13
4.14From this, we get for the determinant of the Hessian
4.15The corresponding expressions for the stable fixed point have to be calculated numerically. The values of *K*^{hyp} and *S*^{hyp} in equation (4.10) require integration of the dynamic equations (4.4)–(4.7) along the heteroclinic orbit connecting the stable and the hyperbolic fixed points. Owing to the nonlinear character of the dynamic equations, these integrations also have to be performed numerically.

To determine the heteroclinic orbit, we used initial conditions such that *M*_{1} and *M*_{2} were located in the stable fixed point of equations (A 7) and (A 8) and the energy was fixed to zero. These conditions left one of the generalized momenta *p*_{1} and *p*_{2} undetermined. Its value was obtained by minimizing the distance of the corresponding trajectory to the hyperbolic fixed point at the origin.

In figure 6, we present the results of our calculation for various cell lengths. The agreement between the simulation results and the approximate expression (4.10) for the switching time is remarkably good taking into account the various approximations we made along the way. The dependence of the switching time *τ* on the total molecule number *Ω* is exponential as indicated by equation (4.10). It also increases exponentially with the cell length *L*. This result is in line with the observation made in [18] that intracellular fluctuations decrease with increasing cell length. The exponential increase in the switching time with the molecule number *Ω* is similar to the behaviour of other biological systems. For example, for genetic toggle switches, an exponential increase in the switching time with the number of transcription factors was found [39,40]. Note, however, that in the case when binding of the transcription factor is not cooperative [41–43] or when transcription of a gene into mRNA is taken into account [44], the switching time depends linearly on protein numbers. Still, in (bio-)chemical systems, exponential dependencies on molecule numbers prevail [21]. To give one more biological example for an exponential dependence, consider the time an evolving population spends in one of different equilibria, which increases exponentially with the (effective) population size [45].

## 5. Discussion

In this work, we have studied the influence of molecular noise on the self-organized polar localization of proteins in rod-shaped bacteria. We found that cooperative attachment to the cell membrane is sufficient to generate asymmetric protein distributions along the bacterial long axis. Using stochastic simulations and an approximate analysis of the weak-noise limit, we found that fluctuations in the protein distribution can lead to switches between two mirror-symmetric solutions. For our analysis, we have assumed protein distributions that are invariant with respect to rotations around the bacterial axis and studied a one-dimensional model. Cooperative binding can spontaneously destroy this invariance, in particular, in short bacteria. It will be interesting to quantitatively investigate the limits of the one-dimensional approach in future work. For large enough ratios between the bacterial length and its diameter, though, we expect the one-dimensional treatment to be appropriate.

The mechanism of self-organization we have analysed is based on cooperative binding of the molecules to a support and their spontaneous unbinding. It is among the simplest schemes generating spatial bi-stability. A possible realization of this system is provided by MinD in *E. coli* in absence of MinE, which is known to produce polar localization [46]. Other molecular mechanisms might be considered [47]. For example, we looked at the case that membrane-bound proteins exist in two states. Proteins bind in a cooperative manner in one state to the membrane and then transit into a second state in which they can unbind from the membrane. In a cellular context, the transition could be linked to the hydrolysis of adenosine triphosphate. If the transition between the two states occurs spontaneously, then the system behaviour is similar to the one discussed above. If proteins in the second state can induce detachment of proteins in the first state, then the switching time will decrease with increasing protein number and system length. This behaviour is similar to the one found for the Min proteins [18,20]. Remarkably, if proteins in the second state only catalyse the transition of proteins in state one to state two, then the system can spontaneously oscillate. Additional cooperative effects can thus have a significant influence on the system behaviour. A detailed account of these systems goes beyond our present aims and will be presented elsewhere.

One might speculate about a possible biological function of the stabilization of polar arrangements with increasing cell length. There are proteins for which a mono-polar localization is essential during cell division, for example, during sporulation of *B. subtilis* and division of *Caulobacter crescentus* [48]. Stochastic switching of these proteins close to or during division would clearly be detrimental. Note, however, that in *B. subtilis* and *C. crescentus* other mechanisms than the one studied in this work are thought to generate mono-polar localization. But then, the Min oscillations become less noisy with increasing cell length, that is, cell age [18]. The reduction of fluctuations might thus be a common feature of bacterial ageing. Clearly, more systems need to be studied from this point of view.

From a theoretical point of view, a number of open questions remain. For example, the effect of cell growth has been studied in the context of polar localization due to nucleoid occlusion [49]. We expect that as long as the rate of cell growth is small compared with the molecular reaction rates and the time proteins diffuse in the cytoplasm before binding to the membrane, then our results should be largely unaffected. Furthermore, we have restricted attention to the weak-noise limit, which is inherent to the use of a Fokker–Planck equation. Clearly, there is a need for techniques to study regimes of stronger noise. For example, one might use the approach developed by Wu & Wang [50] to study the stability of the heterogeneous states as a function of the total protein number and estimate the minimal number of molecules needed to achieve polar localization. In the weak-noise limit, further work is needed to analyse protein structures in three dimensions. In addition to its biological relevance, studying the effects of noise on bacterial protein localization should thus also lead to advances in our understanding of the role of fluctuations in pattern forming systems.

## Funding statement

Funding of this work by Deutsche Forschungsgemeinschaft through the SFB 1027 is gratefully acknowledged.

## Acknowledgements

We thank Vaibhav Wasnik for helpful comments.

## Appendix A. The dynamic equations in the Galerkin approximation

As indicated in the main text, we can stop the expansion in equations (3.19) and (3.20) at *n* = 2,
A1and
A2and still capture the essential features of the system dynamics. Inserting these expressions in equations (3.16) and (3.17) and integrating with respect to *x*, we obtain
A3
A4
A5
A6
A7
A8The coefficients of the diffusion matrix in equation (3.29) then read
A9
A10
A11where *M*_{0}, *C*_{0}, *C*_{1} and *C*_{2} are expressed in terms of *M*_{1} and *M*_{2} by solving equations (A3)–(A6) in steady state.

## Footnotes

One contribution of 13 to a Theme Issue ‘Biophysics of active systems: a themed issue dedicated to the memory of Tom Duke’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.