A one-dimensional three-state run-and-tumble model with a ‘cell cycle’

Abstract We study a one-dimensional three-state run-and-tumble model motivated by the bacterium Caulobacter crescentus which displays a cell cycle between two non-proliferating mobile phases and a proliferating sedentary phase. Our model implements kinetic transitions between the two mobile and one sedentary states described in terms of their number densities, where mobility is allowed with different running speeds in forward and backward direction. We start by analyzing the stationary states of the system and compute the mean and squared-displacements for the distribution of all cells, as well as for the number density of settled cells. The latter displays a surprising super-ballistic scaling \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim t^3$$\end{document}∼t3 at early times. Including repulsive and attractive interactions between the mobile cell populations and the settled cells, we explore the stability of the system and employ numerical methods to study structure formation in the fully nonlinear system. We find traveling waves of bacteria, whose occurrence is quantified in a non-equilibrium state diagram. Grapical abstract


Introduction
Understanding the motion of bacteria has been a classic problem of biophysics [1,2].Bacteria are propelled by their flagellae, whose motor generates a torque which translates into forward or backward motion of the bacteria.The problem has also found interest within the soft matter community, as bacteria are but one example of a much larger class of systems, commonly denoted as microswimmers [3].The run-and-tumble (RT) model of an active particle system is originally motivated by specific features of bacterial motion: this motion only persists for a finite time, the 'run'-time, after which the bacterium stalls, the 'tumble'-period, before continuing its motion typically in a different direction, see e.g.[4].The properties of the basic RT model have been confronted with experiments, e.g. in [5,6].The RT model also relates to other stochastic processes, e.g. the exclusion process [7] or even to the dynamics of quantum particles [8].
RT models in one dimension are a special case within this model class.Here, the bacterium can only switch a e-mail: breoni@hhu.de between left-and right motion in a stochastic manner.One-dimensional RT-models have proven to be an extremely rich field for analytic calculations; exemplary papers dealing with diverse aspects are: confinement [9]; space-dependent velocities, space-dependent transition rates and general drift velocity distributions [10][11][12][13][14]; hard-core particles with spin [15]; inhomogeneous media [16]; attractive/repulsive interactions [17,18]; phase transitions [19]; entropy production [20].Field-theoretic methods have been applied to RT models recently as well [21,22].
In some sense, the (one-dimensional) RT model can be thought of playing in active systems a role analogous to Ising models in equilibrium statistical mechanics.In the very recent past, several works have appeared carrying this analogy further, since they consider the number of 'states' in which the bacterium can find itself to go beyond the dichotomy of left-and right-moving states.Models with three and even more states have been discussed -in our Ising-model analogy, this amounts to looking at active analogues of 'Potts'-type models [23][24][25].
The present paper inserts itself in this line of research by considering a three-state RT model with the states: left-moving, right-moving and sedentary.Our model is motivated by the behavior of the bacterium Caulobacter crescentus (CC), a model organism in microbiology since it has a complex lifestyle [26,27].CC has a bacterial analogue of a cell cycle usually found in eukaryotes; in order to undergo cell division, the bacterium has to switch from its mobile swarmer state to a spatially localized stalked state.Only from the latter state the proliferation of new cells is possible.Our model, capturing this biological feature, is however not limited to CC or bacteria alone.E.g., the green algae Chlamydomonas reinhartii has a similar cell cycle [28] with sedentary and swimming states and also performs a run-and-tumble motion [4].The capacity of cell division in our RT model inks it to the problem of the growth of bacterial colonies.Recently, the authors of [29] developed a growth-expansion model which generates traveling waves in bacterial chemotaxis, in accord with experimental observations.We show that traveling waves also arise in our much simpler 1d run-and-tumble model.
The paper is organized as follows.In Section 2, we introduce our RT-model as a toy model, inspired by the cell cycle of CC.In Section 3, we first focus on the case of free cells for which we derive the conditions for stability of the system when spatial dependencies are neglected.In Section 4 we consider the spatial dependence built into the model and study the mean displacement (MD) and mean squared displacement (MSD) for a single cell in the process of duplicating, both showing a surprising t 3 regime for short times.Allowing the cells to interact via both attraction and repulsion mechanisms, this antagonistic effect is found to lead to structure formation: we numerically find traveling wave solutions of the system density and quantify their occurrence in a non-equilibrium state diagram.Finally we discuss how the model performs with parameter values specific of CC.Section 5 concludes the paper with a discussion of the results of our model and a brief outlook on further work.

The model
Inspired by the reproductive behavior of Caulobacter crescentus we consider a one-dimensional toy model representing bacteria that can actively move rightward, leftward or settle down, and that when settled double in number.We note, that CC performs a run-reverse-flick motion [30], where the bacterium first performs a forward motion, then reverses its direction of motion and in a third step makes a turn mediated by a buckling instability in its flagellum [31].Since our setup is one 1 Graphical representation of the transition rates among different species.These transitions are motivated by the cell cycle of Caulobacter crescentus, that either moves actively or settles down to reproduce.Our model contains three different species: the cells moving to the right ρ + , those moving to the left ρ − and the settled ones ρ 0 .The moving cells can either settle via the rate λ s , move in the opposite direction with λ e or die with µ.Settled cells duplicate via λ d , and generate both a left-and a right-moving cell.
dimensional, the run-reverse-flick motion is equivalent to a run and tumble motion.
The 'cell cycle' of our three-state RT model motivated by CC is summarized in Figure 1.We allow for three populations with the number densities ρ + (x, t), ρ − (x, t) and ρ 0 (x, t), functions of space x and time t, respectively corresponding to right and left movers, and to the sedentary population.The 'cell cycle' step is given by the rate of settling down, λ s , which can occur from either moving state, and the cell doubling with rate λ d with which a sedentary bacterium gives rise to a pair of right-and left-moving cells.The exchange of direction, i.e. the RT step, is denoted by λ e .Finally, µ is the death rate, which we consider for motile cells only.In a proliferating system, this rate prevents exponential growth.This idealized CC-'cell cycle' is implemented in terms of evolution equations for the cell number densities.In the case where there is no death or proliferation, the number densities can also be interpreted as probability densities and the evolution equations correspond to Fokker-Planck equations.
As the bacteria are micron-sized swimmers, we assume a low Reynolds number and overdamped dynamics.To describe this behavior mathematically we first group the three densities into the vector of densities ρ = (ρ + , ρ 0 , ρ − ).The dynamics of the system will then be described by the differential equation which generalizes the standard expression of growthexpansion equations of logistic growth, usually formulated for a single density [29].In Eq. ( 1), the first term is a diffusion term where the matrix D has the form since the sedentary particles do not diffuse.The second term on the right-hand side is a nonlinear diffusion coefficient containing an interaction matrix U of the form The matrix entries describe attractive interactions (negative sign) of the moving cells to regions in which particles have settled and repulsive interactions among settled cells (positive sign) in order to mimic biofilm behaviour.The third term on the right-hand side of Eq.
(1) describes the active motion of the particles in the right and left directions along the line.Hence Finally, we have for the cell cycle or population dynamics part, following the transitions shown in Figure 1, the matrix M given by Given that our run-and-tumble model allows for proliferation and death of cells, it is important to recognize that the population dynamics of Eq. (1) is the linear limit of the more general nonlinear decay-growth equation In Eq. ( 6), M D and M OD are the diagonal and offdiagonal parts of the matrix M, i.e., one has M = M D + M OD .The diagonal part describes the cell number decay, while the off-diagonal part describes the growth of the cell population.In order to limit growth, the non-diagonal term is generally nonlinear and saturating at the carrying capacity, as is common in growthexpansion models, see, e.g.[29].The vector R is thus given by where the carrying capacity is given by the vector The linear limit of Eq. ( 6) is reached for |ρ| |ρ c |.It is important to notice that since R is only applied to one part of the M matrix, the stationary values reached by the population in the linear limit will not necessarily be those given by ρ c .The main benefit of the nonlinear model is that it prevents the number of cells from exploding independently of the parameters.In this manuscript we will mainly focus on the linear case, while explicitly referring to the full nonlinear growth equation if needed.

Free cells
We start by setting the cell interaction parameters κ = κ 0 = 0, and hence consider free cells.

Population dynamics
In this section we further set D = 0 as well as the velocities v + = v − = 0, thus we first study free cells undergoing the pure population dynamics given by This linear system of equations can be solved analytically via matrix calculations, leading to P is the eigenvector matrix of M and E is the diagonal matrix containing the eigenvalues of M, that are We notice that the first two eigenvalues are always negative and therefore stable, while the sign of the third depends on λ s − µ, which can become unstable.This instability facilitates an exponential growth of the colony.In fact, for small values of λ d (λ s −µ) with respect to µ+λ d +λ s the unstable eigenvalue becomes The exponential growth or collapse of the system is therefore decided by the difference of λ s and µ, or in different terms, the separating line between the two behaviors is λ s = µ.It is also worth pointing out that in the case of instant doubling, that is the limit of λ d → ∞, E 3 simply reduces to λ s − µ, as can be seen in Figure 2. Physically this is expected, as in this model cells can double only when settled and can die only when moving, meaning that the growth or decay of the system size depends exclusively on whether a moving cell is faster in settling or dying.In the case of λ s = µ, it is possible Fig. 2 Unstable eigenvalue E 3 (solid lines) as a function of doubling rate λ d for different values of λ s (color code) and µ = 10λ e .The sign of E 3 is the same of λ s − µ, and its value also stabilizes at λ s − µ for very large values of λ d (dashed lines).
to calculate the stationary value of of ρ(x, t → ∞) as a function of the initial conditions ρ(x, 0): ).Since the exchange rate between right ρ + and left ρ − moving cells is symmetric, the amounts of left and right moving cells are the same in the stationary state (ρ + = ρ − , see also Figure 3).Furthermore, if λ d = 2µ = 2λ s all the three populations equilibrate to the same value, independently of the initial conditions.In the case of λ s > µ it is always possible in the frame of the nonlinear growth model to find values of ρ c for which the populations stabilize around the values given by Eq.( 12). Figure 3 shows the linear and nonlinear model equations with different parameters and with the same stationary values.Fig. 3 Space averages of right-ρ + , left-moving ρ − and sedentary ρ 0 cells as functions of time, both for the linear model (solid lines) with λ s = µ and for the nonlinear model (dashed lines) with λ s = 3µ.For both models λ d is set to be equal to λ e , while µ = λ e in the nonlinear model and µ = 2.848λ e in the linear one.As initial conditions we chose the constant values ρ(x, 0) = (0, .1,0.479) for both models.

Density dynamics
We now set the running speeds v ± and the diffusion constant D to finite values, in order to study the evolution of spatial quantities of the system, such as the mean displacement MD = x − x 0 , the mean-squared displacement MSD = (x − x 0 ) 2 and all the higher order moments, where x 0 is the average position of the system at t = 0. Here, the average (•) is defined as ∞ −∞ (•)P (x, t)dx, where the total probability P (x, t) is )dx is the number of cells in phase α and α can be (+, −, 0).
In order to compute averages, we first solve the system by using a Fourier transform (F T ): where ρ(k, t) = F T (ρ(x, t)) is the Fourier transform of ρ(x, t) and k is the wave number conjugate to x.
Similarly to the constant density case, the solution in Fourier space will be given by One can use the solution of this equation to extract the intermediate scattering function (ISF ) The ISF can be related to the different moments of the density [32] by differentiation: valid in one dimension (see Appendix).
We can also define an average for each cell population and the relative ISF : The expression corresponding to Eq.( 17) is then given by First we will discuss the behavior of the whole distribution P (x, t).For simplicity, we will consider the initial condition ρ(x, 0) = (0, δ(x), 0) which is physically relevant, as it describes a cell initially settled in x = 0 in the process of reproducing.We do not focus on the case of an initially mobile cell, as the short-time behaviours of both the MD and MSD turn out to be simply linear and the long-time behaviours are identical to that of the initially settled cell case.We further remark that our analysis does not assume the condition µ = λ s for a stable population in the linear growth model.

Full distribution
When v + = v − , the MD is non-zero and we observe two different regimes: for short times it grows as t 2 , while for long times it is proportional to t, as shown in Figure 4(a).The short-time expansion of the MD in fact yields where The expression shows that both the transition rates and the running speeds have a role in determining this initial scaling regime.This can be interpreted as a composition of the doubling mechanism and the system acceleration given by cells suddenly starting to move.We can further define the typical crossover time t (1) c as the ratio between absolute values of the coefficients of the t 2 and t 3 scalings, as this is the time at which the the t 2 order contribution becomes smaller than the following ones [33,34].This is a good estimate of the average time at which the dynamics is not dominated by the initial doubling anymore: The long-time expansion of the MD yields where Λ is the same of Eq. (10).
As far as the MSD is concerned, in Figure 4(b) we still see a t 2 regime for short times, while the longtime behavior depends on the difference between v − and v + .In the case they are the same, we will only see a diffusive long-time regime while otherwise this diffusive regime transitions into a ballistic one.The smaller the difference between the running speeds, the longer is the time to reach the ballistic regime.We further calculate the short-time expansion of the MSD: where v a = (v 2 + + v 2 − )/2.Again, we define a crossing time t (2) c for the MSD as the ratio between the absolute values of the coefficients of the t 2 and t 3 scalings: If we change the population rates we observe that the growth or decay in the number of cells does not influence qualitatively the scalings we just described for both the MD and MSD.The formula for the long-time expansion of the MSD and the relative crossing time t (2) l between the long-time regimes ∝ t and ∝ t 2 are quite involved, so we refrain from showing them here.Finally, we study directly the full intermediate scattering function F(k, t) as it carries more information than the MSD and MD.In Figure 5 (a), that is in the case of equal velocities, we can see that the real part of F(k, t) that generates the MSD among all other even moments, decays rapidly for small length scales (i.e.large k) while it has three distinct regimes for large length scales.At first the function decays or grows, following the growth in the number of cells, then at time t (2) c it plateaus for a time that grows larger as k gets smaller, and finally decays completely.The plateau, starting after the transition of the cell to its moving stage at time t c , is generated by the active cells going back to the settled stage and not moving anymore, while the final decay represents the long-time diffusive behavior that we have already seen in the MSD.In Figure 5 (b) we see how unequal velocities change the intermediate scattering function by introducing an oscillating behavior at long times.This is a signature of ballistic motion and of a non-vanishing imaginary part of F(k, t) that generates the odd moments like the MD.

Distribution of settled cells
The main feature of the MD and MSD of the settled cells is that they both show an initial t 3 regime, as shown in Figure 6.The short time expansion of the MD is given by with the crossing time between the t 3 and t 4 regimes t c,0 being: The MSD shows the initial t 3 regime as well: with the crossing time t c,0 : The reason why we observe the t 3 -behaviour for short times is the fact that the settled population can only change by doubling, moving and then settling, with each one of these processes being at least of order t.
We also notice that for D = 0 the MSD grows initially with t 4 , as in this case the short time MSD for moving cells grows with t 2 and not t.
The long-time asymptotes for both MD and MSD of the settled particles are identical to those of the whole population.

Attraction to settled regions
We now discuss the case of interacting cells.Our model contains an effective attractive force that pushes the moving cells towards the regions where the density of settled cells is larger.This force is meant to represent how bacteria tend to assemble in resource-rich regions to reproduce or how they accumulate in order to form biofilms [35,36]; therefore the parameter κ > 0 in Eq.( 1).The interaction terms κ∂ x (∂ x (ρ 0 )ρ ± ) render the equation nonlinear such that it is not analytically solvable.Instead we first perform a linear stability analysis around the homogeneous stationary solution to the linear system ρ computed in Eq. ( 12) (see also Fig. 3) by adding a small perturbation δρ(x, t) and neglecting the nonlinear terms in the perturbation (δρ(x, t)) 2 .We then arrive at the following system of equations for the perturbation: where the stationary values for the density are symmetric, ρ+ = ρ− .We apply both a Fourier transform in space and a Laplace transform in time to Eq. ( 30) and solve the resulting characteristic equation of the system.We obtain three different solutions for the eigenvalues of the system s i (k), of which only one, s 1 (k), can have a positive real part.In the following we focus on s 1 (k), since its positive real part introduces instabilities in the system.
First of all, for k → 0, the value of s 1 (k) is one of the eigenvalues of the system matrix where the initial densities are constant, and more specifically the one that can be positive: This means that one of the conditions for the system to be stable is that the number of cells does not grow exponentially, which is expected.
The second limit we consider is k → ∞.We have that This second condition states that the diffusion constant contrasts directly the instabilities generated by a large settling rate and the attractive constant κ, as it disperses too large clusters of active cells, while a large doubling rate helps the stability by reducing the size of groups of settled cells.Knowing the limits of s 1 (k) in (b) c t (2)   l Fig. 4 (a) Mean displacement (MD), (b) mean-squared displacement (MSD), respective crossing times t (1) and shortand long-time approximations for the initial conditions ρ(x, 0) = (0, δ(x), 0)λ e /v + , all rates equal to λ e and D = 0.2v 2 + /λ e .In (b) the solid red line shows unequal swim velocities (v − = 0.9v + ) and the dashed blue line equal swim speeds (v − = v + ).The orange lines represent the short-time approximations, while the green lines are the long-time approximations.(a) 100.0 10.0 1.0 0.1 0.01 t (2)   c t (2)   c t (2)   c t (2)   c t (2)   c 10 2 10 0 10 2 10 4 10 6 c t (2)   c t (2)   c t (2)   c t (2)   c t (2)   l k = 0 and k = ∞, i.e long-and short-range perturbations respectively, we are sure that the system will be unstable if the real part of either of them is larger than zero, giving us two stability conditions for the system: For D = 0, s 1 (k) grows asymptotically like k, making the system always unstable.In Figure 7 we show the behavior of the eigenvalue Re(s 1 (k)) for different values of D. Notice that for the set of parameters considered, if D = 2v 2 + /λ e the stability conditions are only narrowly fulfilled, but the real part of s 1 stays negative for all the values of k.Lastly, when the cell running speeds are not isotropic, the imaginary part of s 1 can be nonzero, meaning that there can be stable periodicity in the system.While the real part of the other two solutions s 2 and s 3 is always negative, their imaginary part is nonzero for large values of k.More specifically, for large k and finite D their imaginary part is proportional to k, while the real part goes with −Dk 2 .A finite imaginary part indicates oscillations in the system, although the negative real part means that these oscillations are only transient.Signatures of these oscillations can also be seen in our numerical solutions (see the next Section).(x(t) x c, 0 t (2)   l Fig. 6 (a) Mean displacement (MD), (b) mean-square-displacement (MSD), respective crossing times t and shortand long-time approximations for settled cells, with initial conditions ρ(x, 0) = (0, δ(x), 0))λ e /v + , all rates equal to λ e and D = 0.2v 2 + /λ e .In (b) the solid red line shows unequal swim velocities (v − = 0.9v + ) and the dashed blue line equal swim speeds (v − = v + ).The orange lines represent the short-time approximations, while the green lines are the long-time approximations.

Repulsion among settled cells
We now include a self-repulsive potential for the cells that do not move, given by κ 0 > 0 in the matrix U in Eq. 1.This repulsion models the need for settled bacteria to not overcrowd any particular region and deplete its resources while reproducing.What is particularly interesting about having both an attractive and a repulsive part in the potential is that the interplay of these two opposing effects can lead to structures forming in the system, as we will show now.If we repeat the analysis described in the last subsection including κ 0 > 0, we find that the limits of s 1 (k) are The main difference with Eqs.( 31), (32) is that s 1 will always be negative for a sufficiently large value of k.This means that if we choose parameters for which s 1 can be positive, its largest root k r will indicate the smallest allowed instability of the system, with size l = 2π/k r .We consequently expect instabilities to form for systems of size L larger than l.As an example of this we numerically calculated the values of k r for different values of the running speeds v + and v − , quantifying their occurrence using two non-dimensional parameters, the maximum speed v m and the reduced difference speed v r defined by We chose specifically to vary the running speeds as they can easily tune the asymmetry of the system, leading to interesting instabilities.In Figure 8 we can see k r as a function of v r and v m , written in units of k 0 = 2π/L.We expect the system to develop instabilities for values of k r > k 0 , so we fitted the separatrix k r = k 0 to a second-order polynomial, v f m (v r ): This particular fit was determined using the linear growth model for the parameter values indicated in the caption to Figure 8.
In order to study the emergence of such instabilities in detail, we further implemented a numerical solver for both Eqs.( 1) and ( 6), using an explicit fourth-order Runge-Kutta algorithm [37] for the time integration and a finite difference scheme in space.We performed calculations both with the linear and the nonlinear growth model, setting respectively λ s = µ and λ s ≥ µ.We use a finite box of length L with periodic boundary and D = 0.001L 2 λ e .In blue we see the parameters for which the system is not large enough to enable instabilities, while in black we have the second order polynomial that fits the k r = k 0 curve.conditions.Setting the time step to ∆t = 10 −4 λ −1 e we calculated ∼ 10 6 steps to ensure that the system settles into a steady state.Our calculations are initialized using the steady-state solutions of the linear system (Eqs.( 12)), to which we add small fluctuations given by Gaussian noise.We find that our system develops wave-like structures, which are static for v + = v − and become traveling waves when v + = v − -see Figure 9 for the linear growth case and Figure 11 for the nonlinear case.Testing different initial conditions, e.g.choosing ρ 0 (x) as a narrow Gaussian peak that approximates an initially settled single cell, we also observed that these wave-like structures always form, even if the specific shape of the wave can be affected.In our analysis we preferred to use the steady-state solution of Eqs.(12) as initial condition, as it makes comparison with the theoretical results of Figure 8 more straightforward.Intuitively, the attractive term κ leads to the formation of peaks, induced by the instability in Eq. (33).These peaks are then stabilized by the repulsive term κ 0 .The asymmetry of the running speeds makes the peaks move.
Migrating bands of bacteria have indeed been observed experimentally [38][39][40][41][42][43] and have also been modeled theoretically [29,44,45], always considering only one species of cells.A particularly surprising feature of our model is that in this final stationary state all three distributions evolve in the same direction at the same speed, independently of the intrinsic running speed of the cells.We replicated the diagram of Figure 8 with numerical integration of the model equation, and the resulting non-equilibrium state diagram is shown in Figure 10.We find a clear transition from a stable system (shown in blue), where all species are constant in space, to the appearance of wave-like structures (shown in red Table 1 Values of the parameters for Caulobacter crescentus taken from [46][47][48][49]. to yellow).The gradient visualizes the change in stationary speed of the waves v s , defined as the speed of the waves in the stationary state divided by √ Dλ e , and is hence non-dimensional.This quantity is almost vanishing near the transition, and grows the further away we move from it.The formation of these waves is typical of systems with a large difference between v + and v − or rather small absolute speeds.We fitted the separatrix to a second order polynomial v f m (v r ) and obtained We find that our numerical calculations and theory are in very good qualitative agreement.

Application to Caulobacter crescentus
Table 1 gives an idea of the experimentally measured values for CC which have been extracted from recent papers on its swimming behaviour [46][47][48][49].It is noteworthy to comment on the running speeds v + , v − .While the torque generated by the flagellar motor differs significantly during forward and backward motion, the resulting velocities are not dramatically different (and, in fact, experimentally hard to measure) [49].We have performed calculations with the parameters of Table 1 for different values of κ and κ 0 which are undetermined from experiments.Since for Caulobacter µ < λ s , we have included the saturating nonlinearity for the growth in the model.The results show that the waves still form provided the ratio κ/κ 0 is large enough (Figure 11).

Conclusions and outlook
In this work we proposed and studied a 1D 3-state model motivated by the cell cycle progression of the bacterium Caulobacter crescentus, including both its run and tumble motion and its reproductive behavior.We first analyzed the free cell space-independent case and calculate the parameter regimes for which the and D = 0.001L 2 λ e .In blue we see the parameters for which the system is stably constant, while in red to yellow we see the parameters for which the system generates traveling wave structures.Examples of both long-time behaviors are shown in their respective area.The gradient shows the stationary velocity of the waves v s , while in black we have the second order polynomial that fits the transition curve v f m .
number of cells grows or declines.Adding the spatial dependence we subsequently determined dynamical quantities of the system such as the mean displacement, the mean-squared displacement and the intermediate scattering function.We found a surprising super-ballistic behavior of the MSD at short times with a t 3 scaling which stems from the interplay of cells doubling and cells starting to swim.
Subsequently, we included attractive and repulsive interactions between cells into our model, representing their tendency to swim towards regions in which cells are settled and to avoid overcrowding.We determined the stability conditions and, using numerical methods, we studied the fully nonlinear system in which we iden-tify traveling waves of cells.Their occurrence is quantified in a non-equilibrium state diagram.
Our model lends itself to further extensions in several ways.E.g., one could account for complex nutrient landscapes and for a more detailed description of the cell cycle, which is well-studied from various aspects [27]; another possible system for application are Chlamydomonas reinhartii cells [28].The cell cycle can be included in cell-resolved simulations such as performed recently in [50,51].Another direction could be a two-dimensional field description that includes the nematic ordering of cells such as in [52].In a higherdimensional model it would also be interesting to see what the effect of different swimming strategies such as run and tumble, run-reverse or run-reverse-flick [30] is.Finally, an exploration of the fully nonlinear model -nonlinear diffusive interactions as well as nonlinear growth -including a full higher-dimensional tumbling behaviour for a multi-species system would be an interesting problem in the context of biofilm growth.1.For the interaction potentials, for which no experimental estimates can be made at present, we chose κ = κ 0 = 10λ −1 e , while for the carrying capacity of the system we set ρ c (x, t) = (0.04, 0.04, 0.04)λ e /v + .
the integral where P (k, t) is the Fourier Transform of P (x, t).Finally, we exchange the order of integration to get Knowing that for the initial conditions that we have chosen ρ(−k, 0) = (0, N (0), 0), we have F(k, t) ≡ P (k, t) P (−k, 0)N (t) = P (k, t)N (t), (41) and hence x n (t) = i n ∂ n P (k, t) ∂k n k=0 = i n N (t)

Fig. 5
Fig. 5 Real part of the intermediate scattering function F (k, t) for (a) equal swimming speeds and (b) unequal swimming speeds (v − = 0.9v + ) for the initial conditions ρ(x, 0) = (0, δ(x), 0))λ e /v + , all rates equal to λ e and D = 0.2v 2 + /λ e .The black lines represent the MSD short crossing time t (2) c and, in the case of different speeds, long crossing time t (2) l .

Fig. 7
Fig. 7 Eigenvalue s 1 (k) as a function of wavenumber k for different values of D, where all rates are equal to λ e , v − = v + and κ = λ −1 e .

Fig. 11
Fig.11Nonlinear growth model, density of left ρ − , right ρ + and sedentary ρ 0 cells as functions of space at different times (increasing from (a) to (c)).Because of the large value of λ e compared to the other rate parameters, the right-moving and left-moving populations have almost the same shape, making the red line disappear under the blue line.The shaded areas indicate the largest peak, and how it moves in time towards the right.We chose as parameters the values typical of CC shown in Table1.For the interaction potentials, for which no experimental estimates can be made at present, we chose κ = κ 0 = 10λ −1 e , while for the carrying capacity of the system we set ρ c (x, t) = (0.04, 0.04, 0.04)λ e /v + .
k) i n ∂ n P (k, t) ∂k n = i n ∂ n P (k, t) ∂k n k=0.