Scale-dependent tipping points of bacterial colonization resistance

Significance The commonly observed stability of natural microbiomes is important for their function, yet the ubiquity of microbiome stability remains enigmatic. The strongest form of stability, colonization resistance, protects residents against invaders and is often associated with specific porous structures, such as skin follicles or intestinal crypts. By systematically probing the colonization of fly gut–derived bacteria in microfluidic pores of varying sizes, we revealed that colonization patterns and invasion rates strongly depend on the pore size. Mathematical modeling shows that bacteria spontaneously tend to organize into a dense colonization-resistant state in pores exceeding a critical size. The scale dependence of stability and resilience could bias ecological filtering in microbiomes and should be considered in the design of microbial ecology experiments.

[2] 49 in terms of the low-density diffusivity D0 and the constant growth rate r. 50 We first expand the cell density as where the frequency ωn of the n th mode is given by 60 ωn = −D0k 2 n + r .
[8] 61 The mode amplitudes an therefore obey an = an(0) exp(ωnt) with the prefactors fixed by the initial conditions.

62
For the empty state to be linearly stable, we require all ωn to be negative, meaning that all mode amplitudes exponentially 63 decay to zero. As the slowest growing mode is n = 0, this implies D0k 2 0 > r, or revealing the establishment transition discussed in the main text. 66 Note that the density decay of non-growing particles (r = 0), which merely diffuse out of the chamber, is on long times 67 controlled by the slowest decaying mode, Φ ∝ exp(−D0k 2 0 t). Since, at the establishment transition, we have D0k 2 0 = r, we see 68 that the diffusive "half-life" of the bacteria just equals their doubling time. This confirms the intuition that the establishment 69 transition occurs when the diffusive outflow is balanced by growth. can lead to negative α (4).
[13] 85 Our discussion of the establishment transition has shown that the frequency of the lowest mode vanishes right at the transition, 86 ω0 = 0, and that the frequencies of all other modes is finite, i.e. ωn = O (1). 87 For ϵ small but finite, we still have that the higher modes have linear relaxation frequencies of order one, ωn = O(1) for 88 n > 1, but the frequency of the lowest mode now assumes small positive frequency of order ϵ, Combining ω0 =O(ϵ), ωn>1 =O(1) with Eq. 12 shows that, at steady state, an>1 with n > 1 is of higher order in ϵ than a0. 92 Thus, to leading order in ϵ, we have Φ(y, t) = a0φ0(y, t) + O(ϵ 2 ), which simplifies the amplitude equation Eq. 12 to 93 ∂tan(t) = ωnan − α 2 D0k 2 n a 2 0 ⟨φn|φ 2 0 ⟩ .

106
To illustrate these results, we consider the example of hard spheres for which α = 1.45 is known exactly (from Eq.

107
(6.12) in Ref. (5)). Our lowest order expansion Eq. 23 then predicts a packing fraction at the floor of the chamber of 108 c(y = 0) = a0 cos(0) = a0 = ϵ6π/(4α) ≈ 3.25ϵ at steady state (t → ∞). Extrapolating from this lowest order expansion, we 109 may estimate that a chamber length of no more than ϵj = 24% above the establishment length is needed for jamming to occur, 110 because then the maximal density at the floor, 3.25ϵj = 0.64, approaches random close packing (Φrcp ≈ 64% for monodisperse 111 spheres). This simple estimate of course ignores non-linear feedbacks, which typically leads to an earlier onset of jamming, as 112 seen in Figs. 2b and 3a.

113
Of particular significance is also the predicted relaxation time to the steady state, the inverse of which is often taken as a 114 measure for the resilience of an ecological system (6). From the exponent in Eq. 23, we see that this relaxation time is given 115 by 1/(2ϵr), which is independent of α and, notably, diverges near the establishment transition. Thus, relaxation can take where we allow for an arbitrary density-dependence in both the collective diffusivity D(Φ) as well as the growth rate b(Φ). As 125 pointed out in the main text, the linear growth rate in our experimental system is to a good approximation constant, which 126 corresponds to b(Φ) = rΦ. The mathematical treatment in this section is independent from this simplifying condition. 127 We define the quantity where we defined an effective potential The boundary conditions imply Π(L) = 0, manifestly ensured through Eq. 26, and ∂yΠ(y = 0) = 0, which we will account for 138 below.

139
Notice that Eq. 28 is formally identical to Newton's equation for a particle at position Π(y) at time y freely falling in a 140 potential U (Π). Hence, we can use the principle of mechanical energy conservation to immediately predict the velocity ∂yΠ of 141 the moving particle in our mechanical analogy, which in our reaction-diffusion problem corresponds to the negative particle 144 manifestly satisfying the no-flux boundary condition ∂yΠ(y = 0) = 0. Note that we use the notation U0 ≡ U (Π0) and 145 Π0 ≡ Π(y = 0).

146
Integrating Eq. 31 over the y yields where  small Φ, we expect D ′ > 0 suggesting that the potential at first rises faster than a parabola. At larger Φ we expect, instead, a 158 negative slope of D(Φ) resulting in a flattening of the potential until jamming leads to rapid rise of the potential.

159
To determine the density at the floor of the chamber, we have to solve the following problem: Let a point mass move down 160 this energy landscape from Φ = Φ0 back down to Φ = 0. The time it takes the moving mass to reach Φ = 0 has to equal the 161 length L of the chamber. If the sojourn time is too small (large) we have to increase (decrease) Φ0. Mathematically, we can 162 formulate this condition using Eq. 34, .
[35] 164 We have thus obtained an equation that can be used to determine a phase diagram as in Fig. 3a from any D(Φ) and b(Φ). 165 Multiple steady states exist if there are multiple Φ0 values with identical sojourn times. The trivial case for which this happens 166 is a simple parabola, which corresponds to the case without density-dependence, D = const. and b(Φ) = rΦ. Then, we have a 167 critical length where any density Φ0 will lead to a marginally stable steady state.

168
More realistically, multiple steady states occur if the potential fluctuates around a parabola. The density interval supporting multiple steady states is bounded by Φ0-values for which [36] This condition can be computed explicitly as follows, where we substituted ξ ≡ Φ0 −Φ in the first line and using Eq. 30 to express the derivative of U in the last line. If Eq. 38 has 169 any solution there must exist multiple steady states. is known to be well described by a Carnahan-Starling equation (9). In terms of the packing fraction, the pressure P (Φ) is Taken together, we obtain the collective diffusion coefficient via where µ0 is the time single particle mobility resulting from the surrounding fluid, η i is a Gaussian random variable with mean discretization. In addition to Brownian motion, the particles move in response to a potential U (r N ) that depends on the full 201 configuration of the system, denoted r N . At any time, there are N total particles, with Nm maturated mother particles and 202 N d growing daughter particles. The interaction potential is decomposable into U r N = U b r N + Ur r N + Uw r N where 203 U b is the bonding potential between mother and daughter particles, Ur is a purely repulsive interparticle interaction and Uw is 204 a confining potential.

205
The bonding potential is taken to taken to be a simple Hookian spring, with stiffness κ, rest length bij, and χij = 1 if particles i and j are a mother daughter pair and χij = 0 otherwise. The specific form of the repulsive interparticle potential is taken to be pairwise decomposable where the pair potential u2(rij), is a WCA potential (10) where ϵ is a characteristic energy scale and where σi is the characteristic size of particle i. All pairs of particles, excluding bonded mother and daughter pairs, interact with these are excluded volume interactions. Finally, the confining potential restricts the the particles motion to approximately an area A = LxLy by imposing a steep potential at x = 0, x = Lx, and y = 0. Specifically, the external wall potential has the form, with a the amplitude of the confining potential and ξ its characteristic length scale.

206
The mechanical system outlined above is conservative and describes the motion of a collection of overdamped particles with a simplified, local description of hydrodynamics (11). At long times, absent added external forces, its evolution would be consistent with thermal distribution and would conserve particle number. In order to model the growth of the bacterial population and division of an individual mother daughter pair, we endow the daughter particles with a time dependent effective size through the deterministic equation of motion with boundary condition σi(0) = 0 and growth rate k. Similarly, to model the budding of the daughter from the mother, the rest distance bij changes in time with the deterministic equation of motion with boundary condition bij(0) = 0 and the same growth rate k. When then impose that when the distance between the 207 mother-daughter pair rij(t) has grown past the effective size of the mother particle, σi, we sever the bond potential by χij → 0 208 between mother-daughter pair ij, provided t > 1/k. An illustration of this criteria is shown in Fig. S12. This growth rule 209 results in an increasing excluded area of the mother-daughter pair that grows approximately linearly in time. All of the results 210 in the main text and below employ this rule with a fixed k. Growth rules that affect in an exponential increase in the excluded 211 area with time yield qualitatively similar results. Further the results presented are for a fixed growth rate k, but generalizations 212 for populations evolving with a distribution of growth rates are also qualitatively similar, provided the distribution is narrow.

213
Synchronized with the bond breaking event, we add new daughter sites to each of the newly divided particles and reinitialize 214 the size and rest length to 0 for each of them. This last step breaks particle number conservation. An illustration of the 215 subsequent increase in Nm is also shown in Fig. S12, which absent mitigating factors will grow exponentially in time.

216
Finally, consistent with the experiments, we apply an absorbing boundary condition at y = Ly. The absorbing boundary 217 condition and particle number growth balance at steady-state, resulting in a mean particle number that depends on the 218 geometry and model parameters. We adopt a unit system where µ = ϵ = k B T = 1, and mature particles sizes σ = 1. This 219 implies that lengths are defined in multiples of the mature particles size, x → x/σ and times in units of the diffusion time for 220 an isolated particle to move its diameter, with D0 = k B T µ0 and accompanying units of time, t → t/(D0σ 2 ). The stiffness κ in 221 principle allows for mechanocoupling between the division and local packing environment, with deterministic division occurring 222 when κ is much larger then the local stress on the mother-daughter pair. In this work, we consider this limit and take κ = 100.

223
The confining potential parameters are taken as a = 10 and ξ = 1.
B.2. Nonequilibrium phase diagram. We have studied the particle based model described above and have found qualitative agreement with both the experimental results and simplified theory presented in the main text. Specifically, we have studied a system 226 with confinement defined by fixed spatial scale Lx = 15 and Ly = 40, seeded with an Nm = 10 initial mother particles. We 227 have found that 10 division cycles is sufficient to relax to a steady state density in the chamber, and evaluated expectation 228 values by averaging over a minimum of 30 additional division cycles. Further, 3-5 independent simulations are used for each 229 expectation value reported. Studies of the effective reaction diffusion model suggests that the function dependence of the 230 system on the length of the chamber enters relative to the establishment length Lest = π D0/4k, so rather then studying 231 different Ly's, we fixed Ly = L and study the dependence on the division rate k, and thus Lest. 232 First, we studied the stead-state density distribution in the chamber. For L/Lest < 1, as expected the density in the chamber is 0 at steady-state on average. For L/Lest > 1, the particles are able to colonize the chamber, and evolve a stationary density distribution. The distribution can be characterized analogously through local packing fraction Φ(y) = ρ(y)π/4 wherê ρ(y) = ⟨δ(y − yi)⟩. Generally, the discrete size of the particles relative to a flat wall will result in an oscillatory density profile when the overall density is larger, ρ > 0.1, which is a result of density correlations induced by their excluded area. Such an oscillatory density profile is not predicted by the simple reaction diffusion model employed in the main text. In order to make contact with that perspective, we report in Fig. S13 density profiles coarse-grained over a the length-scale of the particle(12). We achieve this using by convoluting the number density with a Gaussian, which smooths the profile out. Further, we consider the contribution of the density only from the mother particles, and evaluate 233 expectation values for times that are integer multiples of 1/k.  ⟨[xi(t) − xi(0)] 2 ⟩ < Lx/2, and thus the long time limit strictly goes to 0. We find using Lx = 30 to be large enough that 258 there exists a sufficient separation of timescales between the onset of diffusive motion and its reduction due to confinement that 259 we easily extract a pseudo-time independent self-diffusion constant. This data is reported in Fig. 2. chamber is insulated from the outside by diffusion barrier of descendant cells. Accordingly, the biggest impact of the antibiotic occurs in cavities that, due to the growth rate detriment of the resident, become unjammed and hence invadable (Fig. 4). The 269 larger the growth rate detriment, the larger the range of chamber length that are driven across the tipping point.

270
In the limit of weak selection, and large chamber population size, infiltration happens on a much faster time scale than 271 take over. This allows to analyze the dynamics in two separate steps. The second, takeover, step is familiar from well-mixed 272 populations. The rate of successful take over depends on the competition between selective advantage of the faster growing 273 invader and genetic drift. The first challenge, however, uniquely depends on the spatial structure of the colonized cavities, 274 which is why we mainly focus on the infiltration step.

275
The infiltration step is best analyzed backward in time. As we follow the lineage of a randomly chosen cell backward in time, 276 it is advected towards the floor of the chamber, reflecting the intuitive location advantage discussed above. A balance between 277 self-diffusion and advection leads to a steady state lineage distribution, whose extent scales as the ratio of self-diffusivity and populations. The only difference is that, due to the spatial structure of the population, Ne/N < 1 is mildly and strongly 295 reduced in the gaseous and jammed states, respectively. where Ds is the self-diffusivity and the particle current jp is given by where vp(y, t) is the particle velocity at (y, t) due to any external force.

299
The key part here is the contribution of the lineage current due to self-diffusion. The mathematical form vs = −Ds∂y ln(c) 300 is fixed by the requirement that, for vg = 0, G ∝ c must be a steady sate solution with vanishing current.

301
Suppose, the chamber population has reached a steady state and is large enough so that we can ignore density fluctuations.

302
Assuming there is no external force, vp = 0, the steady state of the ancestor distribution is then given by From this, we can conclude that one particle injected at y, will after relaxation in the cavity take up a fraction This expression is particularly useful for the gaseous phase where gradient-and self-diffusivity are nearly identical. In this case, 309 we have an approximately cosine density profile and, with Dg = Ds, we obtain where c = N/L is the mean population density in the chamber. This expression shows that the neutral invasion success is 312 much less than the well-mixed expectation 1/N because the invading cell has to enter from the supply where the cell density is 313 very low.

314
In the case of a jammed population, the above expressions are not so useful because they require us to know the minute 315 spatial density variations in the jammed phase. Instead, we can exploit the fact that the density hardly varies apart from a 316 boundary layer near the exit of the chamber. A vanishing particle current at steady state requires where we used c(y, t) ≈ c0 in the last step. Hence, the steady state ancestor distribution is a decaying Gaussian, Since Ds/Dg ≪ 1 in the jammed phase, we conclude that fixation becomes small very rapidly when y > λs ≡ 2Ds/r .

321
In our experimental system, we found that λs is about one cell diameter and the thickness of the jammed fraction was 322 about 50 cell diameters, even just after the tipping point. It is, therefore, appropriate to think of the jammed populations as 323 diffusively isolated from the outside environment.   Table S1. Strains cultured in panflute devices. Scale-depedent colonization was observed across multiple species (see Fig. S3), but was not observed for strains which were highly adhesive to walls or had a strong capability of bioflim formation or filamentation. We performed two temperature shift experiments, where we inoculate a panflute at one temperature and, after sufficient relaxation, shift to another temperature, after which we let the system relax again. (Relaxation often took more than 5 hrs.) (a) 30°C to 26°C. Right: The density profiles changed within the different crypts changed substantially the temperature change. Densities are consistently higher at 30°C (orange) than at 26°C (blue). The profiles were measured at steady states with fluorescent microscopy. The insets of the plots show PIV measurements whose slopes indicated the growth rates (see Fig. S4c and d). The growth rate decreased to 87.2 ± 0.8 %. The shift of Lest was analyzed by determined by extrapolating the relation between cavity length and maximal population density at the floor of the crypts to vanishing cell density, similarly to Fig. S5. We found that the establishment length L est increased by 112 ± 11 %. This change was consistent with our theory Lest ∝ 1 √ r

Supplementary Figures
The establishment length shifted upon the temperature change (from 30°C to 22°C). The steady-state density profiles at 30°C (orange) and 22°C (blue) were fitted by a cosine function (black solid and dashed lines, respectively) and normalized. The establishment length was defined by the x-intercept. The relative change of this critical length (6.6 %) was consistent with our theory predicting it to be given by the square root of the relative growth rate change (8.6 %, Fig. S4d). Note that, while these temperature shift experiments are consistent with a pure growth rate change, they come with the caveat that, besides growth rate, additional cell traits might be affected that influence the phase behavior, for instance, the shape of the cells or their intercellular mechanical interactions.