Dividing active and passive particles in nonuniform nutrient environments

To explore the coupling between a growing population of microorganisms such as E. coli and a nonuniform nutrient distribution, we formulate a minimalistic model. It consists of active Brownian particles that divide and grow at a nutrient-dependent rate following the Monod equation. The nutrient concentration obeys a diffusion equation with a consumption term and a point source. In this setting the heterogeneity in the nutrient distribution can be tuned by the diffusion coefficient. In particle-based simulations, we demonstrate that passive and weakly active particles form proliferation-induced clusters when the nutrient is localized, without relying on further mechanisms such as chemotaxis or adhesion. In contrast, strongly active particles disperse in the whole system during their lifetime and no clustering is present. The steady population is unaffected by activity or nonuniform nutrient distribution and only determined by the ratio of nutrient influx and bacterial death. However, the transient dynamics strongly depends on the nutrient distribution and activity. Passive particles in almost uniform nutrient profiles display a strong population overshoot, with clusters forming all over the system. In contrast, when slowly diffusing nutrients remain centred around the source, the bacterial population quickly approaches the steady state due to its strong coupling to the nutrient. Conversely, the population overshoot of highly active particles becomes stronger when the nutrient localisation increases. We successfully map the transient population dynamics onto a uniform model where the effect of the nonuniform nutrient and bacterial distributions are rationalized by two effective areas.


Introduction
Living microorganisms show two distinct types of activity: motility [1][2][3] and proliferation [4][5][6].On the one hand, motile microorganisms locally consume chemical energy, use it to self-propel, and dissipate it due to friction with the fluid surrounding.This flux of energy drives the system out of equilibrium.As a consequence, particles with purely repulsive interaction can show motility-induced phase separation [7][8][9].
On the other hand, proliferation represents a different form of activity.Through cell division and growth, individuals inject biomass into the system, thereby breaking the conservation of mass [4].Mechanical forces due to cell-cell contacts strongly impact the growth process and can be modelled using pressure-dependent growth rates [6,32,33] or a cell cycle impacted by pressure [34,35].Cells more resilient to pressure tend to have a competitive advantage in crowded environments [34,36].The E. coli bacterium, which we have in mind in our model, tends to be less affected by mechanical pressure compared to a variety of epithelial cells [34], so we neglect pressure in this study.However, it should be noted that mechanical self-regulation slows down the elongation rate [37] and in narrow long channels the mechanical forces become the main limiting factor of cell growth [38].Furthermore, for elongated cells the competition between active and passive forces results in a mosaic of locally ordered patches [39].Finally, proliferation and cell death can also liquidise tissue leading to cell diffusion [6] and the transport of cells is additionally amplified during the transition of the colony from two to three dimensions [40].
The interplay between motility and proliferation has already shown interesting dynamics and is crucial in describing and interpreting experimental observations.The combination of density-dependent velocity and growth results in clusters of characteristic size and transient ring patterns [41], a behaviour previously linked to chemotaxis [42][43][44].Furthermore, a system of reaction-diffusion equations modeling growing bacterial colonies in the presence of nutrient [45] demonstrates that motility is crucial to describe the rich variety of spatial patterns observed in experiments [45][46][47].Whereas, immobilizing bacteria reduces variety of possible colony patterns [46].
As soon as cell number is no longer conserved, another aspect becomes relevant -population dynamics.Classic approaches to model it range in complexity from exponential growth for abundant resources, over density-dependent growth rates [5] and coupled differential equations for population-nutrient [48] or predator-prey dynamics [49], to models taking into account stochastic fluctuations [50] or spatial structure [51].The latter has significant impact on population dynamics of fish populations [52], host-parasitoid systems [53], and on the spread of infectious diseases [54].A recent study combined motility and proliferation to model the transition from initial exponential to subexponential growth [55].
In this article, we study the population dynamics of active agents and their emergent collective patterns in a non-uniform nutrient profile using active Brownian particles.Hydrodynamic interactions due to shape-dependent flow fields [56,57] are neglected.
To explore purely proliferation-induced effects, we ignore chemotaxis and its influence on dispersal [58] and clustering [27].For bacteria, this can be realized by genetically suppressing tumbling [59].Our particles or cells divide at a nutrient-dependent rate, while divisions are assumed to be quasi-instantaneous.We consider a constant nutrient influx from a point source.The time scales of motility and proliferation are strongly separated for microorganisms such as E. coli, which also requires large system sizes to explore the non-uniform population and nutrient distributions.Therefore, to keep the computational effort manageable, we squeeze the time scales together.A possible experimental realization is a highly viscous fluid which reduces cell motility [60].
Within such a minimalistic model, the steady population of the bacteria only depends on the ratio between nutrient supply and bacterial death, as expected from a uniform model.In contrast, the nutrient amount is highly governed by its non-uniform distribution around the source, which is tunable by the nutrient diffusion coefficient.In turn, for passive and weakly active particles the nutrient distribution regulates the spatial extent of the bacterial cluster and thereby provides a clustering mechanism that does not rely on chemotaxis or adhesion.However, for high activity these clusters dissolve and the bacteria cover the entire system area even if the nutrient is localized.The transient growth towards a steady state is also heavily influenced by the nutrient diffusion coefficient.We show that it can be mapped onto the dynamics of a uniform model using two limiting cases where we choose effective system sizes for the nutrient and bacterial distribution.
The article is structured as follows.In Sec. 2 we introduce the particle-based and the uniform model.For immobilized bacteria or passive particles, we investigate the steady state in Sec.3.1 and the transient growth dynamics in Sec.3.2.We compare our findings to the results for motile bacteria or active particles in Sec. 4. In the conclusions of Sec. 5 we summarize our results and present an outlook.

Model
In this section, we formulate a particle-based model for dividing bacteria consuming a nutrient, choose feasible bio-inspired parameters, and introduce a model that describes a uniform system as a reference.

Particle-based model
The setup of the two-dimensional system is sketched in Fig. 1 (left).The bacteria are described as overdamped active Brownian particles experiencing thermal noise on their positions and orientations: Figure 1: (left) Sketch of the setup with a source of diffusing nutrients in the centre, reflective boundary conditions, and a quadratic box shape.(right) Mechanical implementation of the cell division.During the division, the circular parent cell is marked in yellow and the daughter cell is marked in blue.Their radii grow in time until the value σ/2 is reached.
with orientation vector u i = (cos φ i , sin φ i ), propulsion velocity v, and translational and rotational diffusion coefficients D T and D R , respectively.For pure thermal noise they are linked to the respective translational and rotational mobilities µ T and µ R by D T /R = µ T /R k B T .For a spherical particle, the Stokes friction coefficients yield D/D R = µ/µ R = σ 2 /3.The Gaussian noises ξ α,i , η i are delta-correlated with zero mean and a variance of one.The volume exclusion force F ij = F σ (r i − r j ) resembles a purely repulsive spring with stiffness k: The Heaviside step function θ(r) cuts off the interaction at the particle diameter σ.
Particles are contained in the system by a reflective boundary condition to avoid accumulation at the wall and an additional repulsive spring interaction between particles and wall to prevent particles from being pushed out of the system area by other particles.We introduce a two-dimensional nutrient concentration field s(r, t).A nutrient is pumped into the system at rate R 0 at location r s , diffuses with diffusion coefficient D N , and is eaten up by the bacteria: The individual uptake rate γg(s(r i , t)) is proportional to the nutrient-dependent growth rate g(s) and the parameter γ is the amount of nutrient required to make a new cell [5].Nutrient amounts will be given in reference to this value.The growth rate follows the Monod function [61] with maximal growth rate g max and characteristic area concentration K s , where the growth rate becomes g max /2.The area concentration is related to the volume concentration, K s = K s σ −1 , used in real three-dimensional systems using the "height" σ of our quasi-two-dimensional system.
To include the nutrient-dependent growth and starvation process, we introduce the rule: In a time interval dt, each particle has the probabilities of dt • g(s(r i , t)) to divide and dt • d to die.Here, we choose a constant death rate d.If a parent cell divides, a daughter cell gets placed on the parent cell.They interact via a repulsive spring force F σ d (t) (r i − r j ), introduced in Eq. ( 2), with growing rest distance σ d (t) = v d t.The speed of the division process is regulated by v d .If the growing rest distance reaches the full particle diameter σ d (t) = σ or the particles drift apart |r i − r j | ≥ σ, the division process is finished.Parent and offspring now interact like all other particles via the force in Eq. ( 2).The division process is sketched in Fig. 1 (right).This implementation is similar to Ref. [35], however, we assume a constant cell size after the cell division is finished.This model is an example of a hybrid particle-continuum reaction-diffusion system.For the simulation we rescale the equations as shown in Appendix A.

Time scales and simulation parameters
A reductionist perspective on our system with relevant time scales is sketched in Fig. 2 (left): I) bacteria and II) nutrient couple via III) growth and IV) nutrient uptake.Additionally, there is V) a constant influx of nutrients and VI) an outflux of bacteria due to the death of individual cells.We define time scales for all six aspects: 3 D R for spherical particles using Stokes friction.

II) nutrient diffusion:
N -time it takes, for example, glucose to diffuse over the characteristic distance.

III) bacterial growth: τ
max -time to consume the nutrient of concentration K s in the local surrounding of the bacteria with area σ 2 .
V) nutrient influx: τ I = γR −1 0 -time to insert the nutrient amount necessary to make a new cell.This can be easily varied in an experimental setup.

VI) bacterial death: τ
The biological system -Escherichia coli in 41 • C water with glucose as the limiting nutrient -dictates the parameters (Appendix B).The resulting time scales span over biology scaling of the exponents scaling of the exponents seven orders of magnitude, as shown in Fig. 2 (right).In our particle-based simulation, this is not feasible.Therefore, we squeeze the time scales together such that the ratios of their orders of magnitude (in units of τ R ), defined as the logarithms, are preserved (e.g.
More concretely, we reduce all the orders of magnitudes by a factor of 3. The procedure is demonstrated for the time scale of bacterial growth: The transformed time scales in Fig. 2 (right) span less than three orders of magnitude and can be realized in the simulation.We recover the simulation parameters from the transformed time scales We additionally vary the propulsion velocity v by choosing a persistence number P = vσ −1 τ R between P = 0 and 20, that corresponds to velocities up to v = 10.562µms −1 .We also implement a sufficiently large spring constant k = 3 × 10 4 k B Tσ −2 .Furthermore, we set the division speed R , which guarantees that the time scale of division σv −1 d is much faster than the time scale of bacterial growth τ G .The simulation box is quadratic with edge length L x = L y = 200 σ.The nutrient source is placed in the centre.The setup is sketched in Fig. 1 (left).At the start of the simulation, N 0 = 0.1•N * bacteria are placed on a square grid and the nutrient field is set to s = 0. Here, N * is the steady-state population in the uniform model introduced in Sec.2.3.The simulation time is 1500 τ R and properties concerning the stationary state use the data between 1000 τ R and 1500 τ R .

Uniform model
As a reference, we introduce a scalar model assuming that nutrients and bacteria are uniformly distributed.When comparing to the simulations, it is of advantage to introduce different areas: the area A S occupied by the nutrients and the area A N ≥ A S of the bacteria.Then, the dynamics of population number N and the total nutrient amount S is described as We assume here a constant influx of nutrient, while outflux of bacteria is only caused by bacterial death.The ratio s = S/A S indicates nutrient density.The ratio A S /A N ≤ 1 comes in since only a fraction of the bacteria, N A S /A N , has access to nutrient and thereby consume them and grow.The case A N = A S presents a limiting case of the population model in Ref. [50] and is distinct from the chemostat, where bacterial and nutrient depletion is caused by material outflux [48].Because the ratio A S /A N can be interpreted as decreasing the growth rate g, we introduce the effective growth rate G(s) = A S A N g(s) = Gmaxs s+Ks with the reduced maximum growth rate G max = A S A N g max .The fixed point of the dynamic equations, gives the steady population N * and nutrient amount S * .Note that N * does not depend on system size.This is plausible because, in steady state the bacterial outflux rate due to death, dN * , and bacterial influx rate due to nutrient-induced division, R 0 /γ, need to balance each other for all system sizes.To characterize the fixed point, we evaluate the eigenvalues of the Jacobian, where the parameter relates the nutrient area A S to the critical area A c N * .The fixed point is always stable because Re(λ 1,2 ) < 0, but the critical nutrient area A S = A c N * corresponding to Ω = 1 separates overdamped (Ω < 1) and damped oscillatory (Ω > 1) population dynamics from each other.In Fig. 3 we illustrate the dependence of Ω on A S and A N .The red line indicates Ω = 1.While for A S = A N the dynamics switches between overdamped and damped oscillatory, the parameter Ω is always larger than 1 to the right of the red line, for example, for A N = A box .
In the following, two cases will be relevant: (1) nutrient and bacteria cover the same area A N = A S and (2) bacteria cover the entire system area A N = A box while nutrient can be localized with A S ≤ A box .We will use case 1 to fit the simulation data of passive particles, while case 2 applies to highly active particles that always spread over the whole system due to their activity.Figure 4 illustrates that the population dynamics of the two cases show a different trend when the area of nutrient A S is increased.For A N = A S (left) an increase of A S increases the characteristic overshoot in N , because the lower densities result in a weaker coupling between nutrient and bacterial population.For A S ≤ A N = A box (right), the characteristic overshoot in N always exists, but now decreases with an increase of A S because a larger fraction of the bacteria can interact with the nutrient and coupling is stronger.Furthermore, with decreasing A S the dynamics towards steady state is considerably slowed down.

Results for Passive Particles
We start our investigations with passive particles (P = 0), which helps us to work out the effect of activity or self-propulsion in Sec. 4. Figure 5 gives an overview of our essential results from the 2D simulations.For a small nutrient diffusion coefficient (left), the population N grows continuously until the steady state is reached, while the nutrient amount S goes through a maximum.The population strongly clusters around the source of the nutrient, which is nonuniformly distributed.In contrast, for large D N = 10 000 σ 2 τ −1 R (right), the nutrient spreads evenly over the whole system area and the population is almost uniform during the growth process.The population relaxes faster towards the steady state for the small D N , but the steady population is not affected by the nonuniform distribution of nutrient and particles.
In the following, we will establish a connection between the population dynamics and heterogeneities in the system.We will introduce an effective system size within the uniform model of Eq. ( 6) to understand the role of such heterogeneities.Following the observation that nutrient and bacteria cover roughly the same area, we equate both areas in the uniform model, A N = A S .This will give us insights into the steady state (Sec.3.1) and the transient phase (Sec.3.2).The parameters not explicitly mentioned are chosen as discussed in Sec.2.2.

Steady State
We start with analysing the steady state.Figure 6 shows snapshots of the particle distribution with different degree of heterogeneity tuned by the nutrient diffusion coefficient D N .For small D N the nutrient distribution is strongly peaked and overlaid by the strongly clustered particles that only occupy a fraction of the system area.Increasing D N allows the population to spread over larger areas.For D N = 10 000 σ 2 τ −1 R the cluster has completely dissolved and the entire area is covered by particles.In other words there is an effective system area A eff occupied by the particles and nutrient that grows with the nutrient diffusion coefficient and reaches the total system area for sufficiently large D N .First, we look at population N * sim and nutrient amount S * sim in the steady state, which we determined in the 2D simulations.Figure 7 shows both quantities plotted versus D N for different nutrient input rates R 0 .The steady state of the uniform model in Eq. ( 7) with N * = R 0 /(dγ), S * = A S G −1 (d), and A S = A box is represented by the dashed lines.Most notably, the population in the steady state, N * sim , is independent of D N and, thus, it is not influenced by the degree of heterogeneity, for example, in the nutrient distribution.The argument used in Sec.2.3 explaining the independence of N * from system size also applies here.Bacterial influx R 0 /γ and bacterial outflux dN * always need to balance for all heterogeneities and system sizes.In contrast, the nutrient amount strongly depends on D N and only for D N /σ 2 τ −1 R → ∞ converges to the prediction of the uniform model, for which we used the system size, A S = A box .This is the case, because for large D N the nutrient is distributed almost uniformly and the uniform model applies with Since S * ∝ A S , one could be tempted to introduce an effective system size A eff to explain the results from the particle-based simulation in Fig. 7, right, which deviate from the prediction of the uniform model for small and intermediate D N .However, we already know that the occupied area A eff grows with D N and, therefore, would expect S * sim to grow with D N .But this is only the case for D N larger than ca. 10 2 σ 2 τ −1 R , for smaller D N the nutrient amount in the steady state rises again.
We propose that this is caused by the limited uptake rate per area the particles can provide.We consider the specific case R 0 = 3000dγ = 90 τ −1 R γ.The diffusion equation Eq. ( 3) of the nutrient is solved on a lattice with constant 2σ.Thus, the nutrient from the delta-peaked source is initially distributed over the area 4σ 2 , which results in an input rate per area of 22.5 τ −1 R σ −2 γ.In contrast, the maximum uptake rate per individual, γg max , and the maximum number density in hexagonal tight packing, [62], limits the maximum uptake rate per area to  7) and the full lines in the right plot are guides to the eye.The inset compares steady nutrient amount versus D N for particles with and without hard-core repulsion for R 0 = 3000dγ.
This value is much lower than the input rate per area at the position of the source.The nutrient needs to diffuse away before it can entirely be consumed, which needs more time for smaller D N .As a consequence, the nutrient amount in the steady state rises.To validate this hypothesis, we performed simulations with particles, which do not interact by volume exclusion.Thus, the maximum density and thereby the maximum uptake rate per area is no longer limited.Indeed, the inset in Fig. 7, right shows that now the steady nutrient amount S * sim is monotonously rising with D N , as predicted.

Transient dynamics
Now that we established a first connection between nutrient diffusion coefficient D N , the nonuniform steady-state nutrient distribution, and the area occupied by the particles, we proceed by looking at the transient dynamics.First, we look closer at the spatial distribution of the growing and evolving particle population.The snapshots in Fig. 6 convey that in steady state the particle distribution is radially symmetric.This is certainly not the case in the beginning of the growth process as the snapshot (a) in Fig. 5 already shows.In Fig. 8 we consider the system dynamics for N * = 3000 and D N = 200 σ 2 τ R .However, to amplify the effect of statistical fluctuations, a smaller initial population N 0 = 100 is chosen in contrast to N 0 = N * /10 used otherwise.The simulation starts with a small population (a).The nutrient quickly builds up until it reaches a maximum (b); now particle growth is possible in large parts of the system.Individual particles establish clusters across the entire system, the most prominent around the nutrient source (c).Because stochastic fluctuations strongly affect the distribution of the few initial particles, the clusters are unevenly distributed and also the nutrient distribution becomes slightly asymmetric.The clusters grow, consume nutrients, and overshoot the carrying capacity of the steady state (d).In parallel, a depletion zone of nutrients around the central cluster develops (c).As a consequence of the nutrient depletion, the clusters in the outer regions start to dissolve (e).As the system relaxes further toward the steady state, the outer clusters completely dissolve, the central cluster becomes symmetric, and the depletion zone vanishes (e,f).
For larger D N the asymmetry tends to be more prominent, while for small D N only the centre is covered by nutrients and clusters cannot grow in the outer regions.However, the nutrient depletion zone is more pronounced for small D N , as Fig. 5 clearly shows.The diffusion coefficient D N = 200 σ 2 τ R depicted in Fig. 8 is an intermediate value, where both effects can be observed.
We further rationalize the relation between population dynamics and nonuniform nutrient distribution.In Fig. 9 (left) we plot population N versus time for different D N .The population dynamics is more strongly damped for more heterogeneous nutrient distributions, i.e., when D N is small.However, for increasing D N we see that a clear overshoot in N develops.We did already observe the same transition from overdamped relaxation to damped oscillations in the uniform model for increasing and equal areas covered by nutrient and bacteria, A S = A N , as shown in Fig. 4 (left).In the following we will introduce an effective system size A eff = A S = A N to model this transition in the particle-based simulations and thereby account for heterogeneities in the system.
The grey dashed lines in Fig. 9 (left) show fits of the simulation results using the uniform model of Eq. ( 6).To fit the curves, we use the fit parameter A eff = A N = A S that can deviate from the simulation box size A box .For D N /σ 2 τ −1 R = 10, 1000, 10000 the fit captures the dynamics well, only for D N = 100 σ 2 τ −1 R a clear deviation is apparent.So, by choosing an effective system size A eff that accounts for the fact that particles and nutrient display a nonuniform distribution in the simulation box, we are able to model the simulated dynamics.The fitted effective system area A eff over the diffusion coefficient is shown in Fig. 9 (right) for different steady populations N * .We observe a monotonous relation between heterogeneity (D N ) and effective system size.In particular, for smaller N * the effective system size jumps at around D N = 100σ 2 τ −1 R , while for large D N it approaches the total system area.This all agrees well with our former observation that for small D N the particle cluster only occupies a fraction of the total system, while for large D N it covers the entire area.

Results for Active Particles
So far we looked at immobilized bacteria modelled as passive Brownian particles.In this part we demonstrate how motility changes both steady-state distribution and population dynamics.We model motile bacteria as active Brownian particles using Eq.(1).First, the case of high activity with propulsion velocity v = 20 στ −1 R corresponding to a persistence number P = 20 is explored.We can still describe its population dynamics with the uniform model with modified choices for the areas A N and A S .Second, we consider a lower activity with v = 10 στ −1 R and P = 10.Here, the population dynamics is intermediate between the highly active and the passive case.

High Activity P = 20
Figure 10 (left) shows the motion of a particle (blue) and all its descendants.The spacetime trajectories known for conserved particle number are replaced by space-time trees of lineages of particles being born ( * ) and dying (+) as discussed in Ref. [4].The lineage depicted starts with one particle (blue trajectory) having two direct offspring (orange, green) with one of them having another offspring (red).Then, the lineage goes extinct.We note when nutrient is localized around the centre, the birth of active particles is also limited to the central region, but their motility allows the particles to explore the area and die more evenly distributed over the system.This is illustrated in Fig. 10 (centre, right).This feature is also apparent in Fig. 11, where we show the snapshots of the steady state for different diffusion coefficients D N for N * = 3000.In strong contrast to the proliferation-induced clusters observed for passive particles in Fig. 6, the distribution of highly active particles is almost uniform for all D N , with slightly higher particle densities around the source for low D N .For small D N the nutrient is localized, but the peak is less pronounced than for passive particles.This can be explained by the absence of the proliferation-induced particle cluster around the nutrient, which decreases the immediate uptake and allows the nutrient to diffuse further away from the source.
The population dynamics of these highly active particles is presented in Fig. 12 and shows a different trend than for passive particles, which we addressed in Fig. 9 (left).While for passive particles the population overshoot increases with D N , for active particles the overshoot is more pronounced for small D N .The same two trends are observed in the scalar population model.We obtained them in Fig. 4

by setting either A
How can we transfer this insight to the particle-based model?The assumption made for passive particles that bacteria and nutrient cover roughly the same areas, A N = A S , does not hold for active particles.Instead, we use the observation that bacteria spread approximately evenly over the entire system area and set A N = A box , while we let the area of the nutrient A S vary with the diffusion coefficient D N .Now, if nutrient is more localized, a smaller fraction of the bacteria interacts with the nutrient.Thus, the coupling between both species is weaker and the overshoot in the particle population is stronger.
To further validate this scenario for active particles, we fit the population model to the population curves obtained from the simulations using the area of the nutrient A S = A eff as the fit parameter.The resulting curves, gray dashed lines in Fig. 12, show good agreement with the particle-based simulation data.In particular, the time positions of the overshoots are always well captured.The fitted area of the nutrient A eff grows with the diffusion coefficient and approaches A box for larger D N as we expected.Only for low diffusion coefficients D N = 2 and 5 σ 2 τ −1 R significant deviations are observed around the overshoots.Here, transient clusters occur and therefore the assumption that particles are evenly distributed does no longer apply.The forming and dissolving clusters also result in a time-dependent effective size A N of the bacterial distribution, which is not included in the model.Nonetheless, the tendency of a stronger overshoot for smaller D N , the time positions of the overshoots, and the monotonous relationship between D N and A S are all well captured.

Low Activity P = 10
So far we discussed passive and strongly active particles with distinct collective behaviour and population dynamics.Now the intermediate case of low activity (persistence number P = 10) is illuminated.Snapshots of the steady state in Fig. 13 (top) show that for small D N clusters form, but are less pronounced compared to the passive case in Fig. 6, and a significant fraction of particles is outside the central cluster.Clusters dissolve for intermediate diffusion coefficients, e.g., D N = 100 σ 2 τ −1 R but the bacterial density is still higher around the source.
The transient population dynamics is shown in Fig. 13 (bottom).Similar to the highly active case, the population overshoot is more pronounced for small D N .Fitting the observed dynamics with the uniform model and A N = A box as before (gray dashed line), works very well for the largest D N = 1000 σ 2 τ −1 R .However, already at D N = 100 σ 2 τ −1 R and very pronouncedly at D N = 20 σ 2 τ −1 R the fitted area for the nutrient, A S > A box , is unphysical.The reason is that bacteria are localized such that the premise of evenly distributed bacteria in the whole system, A N = A box , which we used for highly active particles, does not apply.As an alternative, we treated A S and A N as independent parameters.However, this generated ambiguous parameter pairs.In the end we tried to fit with A S = A N (black dashed lines), which gives realistic values A S = A N < A box and a very good fit for D N = 100 σ 2 τ −1 R .The fit clearly does not reproduce the height of the overshoot for D N = 5 σ 2 τ −1 R since the assumption of a uniform distribution of the active particles is strongly violated.

Conclusion
We presented a minimalistic model for bacterial growth in a nonuniform nutrient environment, where we couple proliferating active or passive Brownian particles to a diffusing nutrient emitted from a point source.This allows us to gain insights into the interplay between collective behaviour and population dynamics of bacteria.The heterogeneity of the nutrient distribution around the point source highly depends on its diffusion coefficient D N .To keep the computational effort manageable, we needed to squeeze the time scales of motility and growth together.
Using numerical simulations, we show that localizing the nutrient supply limits the effective area of bacterial growth and results in proliferation-induced bacterial clustering around the nutrient source for passive and weakly active particles.This presents an alternative clustering mechanism that does not rely on chemotaxis or adhesion.In contrast, highly active particles can disperse over the whole system during their lifetime resulting in a nearly uniform distribution.Interestingly, the space-time trajectories of non-dividing particles are replaced by space-time trees of lineages that can persist, flourish, and die out [4].For both active and passive particles the steady population is purely determined by the ratio of nutrient influx to cell death and can be predicted by a uniform model.
The transient dynamics of the system strongly depends on the nutrient and bacterial distribution, which are determined by the nutrient diffusion coefficient and bacterial motility or activity.For passive particles and rather uniform nutrient profiles at larger D N , fluctuations in the small initial bacterial distribution amplify in the transient phase because cells can form clusters throughout the entire system.This results in an asymmetric bacterial distribution during the transient phase, which dissolves when approaching the steady state.In a peaked nutrient distribution occurring at small D N , bacterial growth is limited to a small region during the entire transient phase and the bacterial distribution looks rather symmetric around the nutrient source.
Interestingly, in non-motile systems the bacterial swarm and the nutrient roughly occupy the same area, which influences population dynamics.If the nutrient and thereby the bacterial cluster are strongly localized, the population number quickly approaches the steady state.We propose this is due to the localized nutrient almost entirely covered by bacteria, which results in a strong coupling.At increased nutrient diffusion, where the profile is more uniform, the nutrient strongly overshoots and is followed by a population overshoot.Here, nutrient and bacteria cover a larger area and the coupling is weaker.
In contrast, highly active particles cover the entire system area even for peaked nutrient distributions and the trend reverses.For strongly localized nutrient at small diffusion constant the bacterial population strongly overshoots because a large fraction of the bacteria do not interact with the nutrient, which means weak coupling.If the nutrient profile becomes flatter, the coupling is stronger since a larger fraction of the bacteria can interact with the nutrient.
The observed difference in the dynamics can be rationalized by a uniform population model, where we introduce separate areas of the bacterial (A N ) and nutrient (A S ) distributions.Using A N = A S for the passive case and A N = A box for the case of high activity, the population dynamics can be well fitted with a single parameter A S .The relation between diffusion coefficient and effective nutrient area is monotonous in both cases.In other words, a more heterogeneous nutrient distribution results in a smaller effective system size.For the intermediate case of weakly active particles, the uniform model has limited applicability in the case of slow nutrient diffusion, where both limiting cases of zero and high activity cannot be applied.
We demonstrated that our minimalistic model allows us to build a bridge between analytic population models and particle-based simulations and to highlight differences between dividing active and passive particles.The model makes general statements that serve as a guideline for future theoretical and experimental investigations.Our insights into the role of nutrient heterogeneity can help to understand how an environment influences the individual reproductive success [63][64][65][66], the spread of neutral mutations [67,68], and the population dynamics in changing environments [69].To make quantitative predictions, our model can be extended to include processes such as chemotaxis [27][28][29][30][31], density/pressure-dependent growth rates [6,32,33], cell elongation [37,39], polydispersity [70,71], and more complex setup geometries.

Figure 3 :
Figure 3: Parameter Ω over nutrient area A S and bactrial area A N ≥ A S .For G max < d death dominates growth and the population dies out.Parameters are chosen as discussed in Sec.2.2 with steady population N * = 3000.

Figure 4 :
Figure 4: Population dynamics for (left) case 1 with A N = A S and (right) case 2 with A S < A N = A box for varying A S .Areas are given in units of the system area A box .In both cases a characteristic overshoot of the bacterial population is observed.Parameters are chosen as discussed in Sec.2.2 with steady population N * = 3000.

Figure 5 :
Figure 5: (Bottom) Dynamics of population N and nutrient amount S with nutrient input rate R 0 = 3000dγ for the diffusion coefficients D N = 10 σ 2 τ −1 R (left) and D N = 10 000 σ 2 τ −1 R (right).(Top) The snapshots show different stages of the growth process as indicated in the bottom plots.Dots represent particles and the green shading the nutrient field.For small D N the nutrient is almost entirely covered by particles in the steady state.

Figure 7 :
Figure 7: (left) Steady population of the particle-based model, N * sim , and (right) steady nutrient amount S * sim (in units of characteristic nutrient amount γ) plotted versus nutrient diffusion coefficient D N for different nutrient input rates R 0 (in units of dγ).Dashed lines indicate the predictions of the uniform model, N * and S * , as stated in Eqs.(7) and the full lines in the right plot are guides to the eye.The inset compares steady nutrient amount versus D N for particles with and without hard-core repulsion for R 0 = 3000dγ.

Figure 8 :
Figure 8: Transient dynamics of a system with steady population N * = 3000, initial population N 0 = 100, and diffusion coefficient D N = 200 σ 2 τ R .Left: Snapshots of evolving particle population and nutrient distribution (green shading) at different times as indicated in the right plot.The green curve illustrates the nutrient concentration along the grey dashed line.Right: Population N and nutrient amount S plotted versus time.

Figure 9 :
Figure 9: (left) Population dynamics of the particle-based simulations (green) with N * = 3000 for different diffusion coefficients of the nutrient, D N .The dynamics is fitted with the uniform model of Eq. (6) assuming A N = A S = A eff (grey) using A eff as a fit parameter.(right) Relation between the diffusion coefficient of nutrient D N and the effective system area A eff for different steady populations N * .

Figure 10 :
Figure 10: (right) Space-time tree of a lineage starting with one active particle.The stars and crosses indicate cell division and death, respectively.(centre) Spatial probability distribution of birth/cell divisions and (right) particle death.The data is obtained for active particles with persistence number P = 20, steady population N * = 3000, and a small nutrient diffusion coefficient D N = 10 σ 2 τ −1 R .

Figure 13 :
Figure 13: Top: Snapshots for active particles with lower activity (P = 10) for different D N with N * = 3000.Bottom: Population dynamics.The gray dashed lines show fits with the uniform model with A N = A box as for highly active particles, while the black dashed lines use A N = A S as in the passive case.