Modelling microbiome recovery after antibiotics using a stability landscape framework

Treatment with antibiotics is one of the most extreme perturbations to the human microbiome. Even standard courses of antibiotics dramatically reduce the microbiome’s diversity and can cause transitions to dysbiotic states. Conceptually, this is often described as a ‘stability landscape’: the microbiome sits in a landscape with multiple stable equilibria, and sufficiently strong perturbations can shift the microbiome from its normal equilibrium to another state. However, this picture is only qualitative and has not been incorporated in previous mathematical models of the effects of antibiotics. Here, we outline a simple quantitative model based on the stability landscape concept and demonstrate its success on real data. Our analytical impulse-response model has minimal assumptions with three parameters. We fit this model in a Bayesian framework to data from a previous study of the year-long effects of short courses of four common antibiotics on the gut and oral microbiomes, allowing us to compare parameters between antibiotics and microbiomes, and further validate our model using data from another study looking at the impact of a combination of last-resort antibiotics on the gut microbiome. Using Bayesian model selection we find support for a long-term transition to an alternative microbiome state after courses of certain antibiotics in both the gut and oral microbiomes. Quantitative stability landscape frameworks are an exciting avenue for future microbiome modelling.

as a 'stability landscape': the microbiome sits in a landscape with multiple stable equilibria, 23 and sufficiently strong perturbations can shift the microbiome from its normal equilibrium to 24 another state. However, this picture is only qualitative and has not been incorporated in 25 previous mathematical models of the effects of antibiotics. Here, we outline a simple 26 quantitative model based on the stability landscape concept and demonstrate its success 27 on real data. Our analytical impulse-response model has minimal assumptions with three 28 parameters. We fit this model in a Bayesian framework to previously published data on the 29 year-long effects of four common antibiotics (ciprofloxacin, clindamycin, minocycline, and 30 amoxicillin) on the gut and oral microbiomes, allowing us to compare parameters between 31 antibiotics and microbiomes. Furthermore, using Bayesian model selection we find support 32 Introduction 40

Stability and perturbation in the microbiome 41
The human microbiome is a complex ecosystem. While stability is the norm in the gut 42 microbiome, disturbances and their consequences are important when considering the 43 impact of the gut microbiome on human health (1). A course of antibiotics is a major 44 perturbation, typically leading to a marked reduction in species diversity before 45 subsequent recovery (2). Aside from concerns about the development of antibiotic 46 resistance, even a brief course can result in long-term effects on community composition, 47 with species diversity remaining lower than its baseline value for up to a year afterwards 48 (3). However, modelling the recovery of the microbiome is challenging, due to the difficulty 49 of quantifying the in vivo effects of antibiotics on the hundreds of co-occurring species that 50 make up typical microbial communities within the human body. 51 Artificial perturbation experiments are widely used to explore the underlying dynamics of 52 macro-ecological systems (4). In the context of the gut microbiome, the effects of 53 antibiotics have previously been investigated descriptively (5-7). However, despite interest 54 in the application of ecological theory to the gut microbiome (8), there has been limited 55 quantitative or mechanistic modelling of this response. In general, the diversity of the 56 microbiome falls before recovering, but the nature of this recovery remains unclear. While 57 responses can appear highly individualized (7) this does not preclude the possibility of 58 generalized models applicable at the population level. 59 Applying mathematical models to other ecological systems subject to perturbation has 60 given useful insight into their behaviour (9-11). Crucially, it allows the comparison of 61 different hypotheses about the behaviour of the system using model selection. Developing In one popular schematic picture taken from classical ecology, the state of the gut 109 microbiome is represented by a ball sitting in a stability landscape (1,(23)(24)(25). Perturbations 110 can be thought of as forces acting on the ball to displace it from its equilibrium position 111 (25) or as alterations of the stability landscape (26). While this image is usually provided 112 only as a conceptual model to aid thinking about the complexity of the ecosystem, we use 113 it here to derive a mathematical model. 114 We model the effect of a brief course of antibiotics on the microbial community's 115 phylogenetic diversity as the impulse response of an overdamped harmonic oscillator 116 ( Figure 1; see Methods), and compare parameters for four widely-used antibiotics by fitting 117 to empirical data previously published by Zaura et al. (3). This model is significantly less 118 complicated than other previous models developed for similar purposes, but still captures 119 some of the essential emergent features of such a system while avoiding the 120 computational difficulties of fitting hundreds of parameters to a sparse dataset. After 121 demonstrating the effectiveness of this modelling approach for the gut and oral 122 microbiomes, we also show that the framework can easily be used to test hypotheses 123 about microbiome states. We compare a model variant which allows a transition to a new 124 equilibrium, and find that this model is better supported for clindamycin and ciprofloxacin, 125 allowing us to conclude that these antibiotics can produce state transitions across different 126 microbiomes. This modelling approach can be easily applied to sparse datasets from 127 different human microbiomes and antibiotics, providing a simple but consistent 128 foundational framework for quantifying the in vivo impacts of antibiotics. 129 Taking inspiration from classical ecological theory, the microbiome can be considered as 132 an ecosystem existing in a stability landscape: it typically rests at some equilibrium, but 133 can be displaced ( Figure 1A). Any quantitative model of the microbiome based on this 134 concept requires a definition of equilibrium and displacement. While earlier studies sought 135 to identify a equilibrium core set of 'healthy' microbes, disturbances of which would 136 quantify displacement, it has become apparent that this is not a practical definition due to 137 high inter-individual variability in taxonomic composition (25). More recent concepts of a 138 healthy 'functional core' appear more promising, but characterization is challenging, 139 particularly as many gut microbiome studies use 16S rRNA gene sequencing rather than 140 whole-genome shotgun sequencing. For these reasons, we choose a metric that offers a 141 proxy for the general functional potential of the gut microbiome: phylogenetic diversity (25). 142 Diversity is commonly used as a summary statistic in microbiome analyses and higher 143 diversity has previously been associated with health (27) and temporal stability (28). Of 144 course, describing the microbiome using only a single number loses a great deal of 145 information. However, if we are seeking to build a general model of microbiome recovery 146 after perturbation, it seems appropriate to consider a simple metric first to see how such a 147 model performs before developing more complicated definitions of equilibrium, which may 148 generalise poorly across different niches and individuals. 149 We assume the equilibrium position to have higher diversity than the points immediately 150 surrounding it (i.e. creating a potential well) ( Figure 1B). However, there may be alternative 151 stable states ( Figure 1B) which perturbations may move the microbiome into ( Figure 1C). 152 These states may be either higher or lower in diversity; for our purposes, all we assume is 153 that they are separated from the initial equilibrium by a potential barrier of diversity i.e. a decrease of diversity is required to access them, which is what keeps the microbiome at 155 equilibrium under normal conditions. 156

The model 157
Mathematically, small displacements of a mass from an equilibrium point can be 158 approximated as a simple harmonic oscillator (29) for any potential function (continuous 159 and differentiable). This approximation comes naturally from the first terms in the Taylor 160 expansion of a function (30), and can be extremely accurate for small perturbations. By 161 assuming the local stability landscape of the microbiome can be reasonably approximated 162 as a harmonic potential, we are assuming a 'restoring' force proportional to the 163 displacement x from the equilibrium position (െ݇‫)ݔ‬ and also a 'frictional' force acting 164 against the direction of motion (െܾ‫ݔ‬ሶ ). The system is a damped harmonic oscillator with the 165 following equation of motion: 166 Additional forces acting on the system -perturbations -will appear on the right-hand 168 baseline and fitting the model (taking into account the whole continuous temporal 207 response rather than pairwise comparisons of absolute diversity) it appears that slow 208 reconstitution cannot be the whole story. Instead, the skewed distribution of residuals after 209 a year, when the response has flattened off, indicates that the longer-term dynamics of the 210 system do not obey the same impulse response as the short-term dynamics. A scenario 211 involving a long-term transition to an alternative stable state is consistent with this 212 observation ( Figure 1). We therefore developed a variant of the model to take into account 213 alternative equilibria, aiming to test the hypothesis that the microbiome had transitioned to 214 an alternative stable state. 215

A model allowing antibiotic-induced state transitions 216
In our approach, a transition to an alternative stable state means that the value of diversity 217 displacement from the original equilibrium will asymptotically tend to a non-zero value. 218 There are many options for representing this mathematically; for reasons of model There were other clear differences in response between antibiotics. For example, 263 clindamycin caused a decrease in the anaerobic Gram-positives Ruminococcaceae after a 264 month, whereas ciprofloxacin had no effect. There was also an individualized response: 265 ciprofloxacin led to dramatic increases in Erysipelotrichaceae for some participants, and for these individuals the increases coincided with marked decreases in Bacteroidaceae, 267 suggesting the relevance of inter-family microbial interactions (Supplementary File 1). 268 Comparing relative abundances at the family level, there were few differences between 269 community states of different treatment groups after a year. Equal phylogenetic diversity 270 can be produced by different community composition, and this suggests against consistent 271 trends in the long-term dysbiosis associated with each antibiotic. However, we did find that 272 Peptostreptococcaceae, a member of the order Clostridiales, was significantly more 273 abundant in the clindamycin group when compared to both the ciprofloxacin group and the 274 placebo group separately (p < 0.05, Wilcoxon rank sum test). In a clinical setting, 275 clindamycin is well-established to lead to an increased risk of a life-threatening infection 276 caused by another member of Clostridiales: Clostridium difficile (34). Long-term reductions 277 in diversity may similarly increase the risk of overgrowth of pathogenic species. 278

Connection to generalized Lotka-Volterra models 279
We sought to establish a link between our framework and the conventional 'bottom-up' 280 approach of generalized Lokta-Volterra models. We investigated the behaviour of a 3 281 species Lotka-Volterra system to establish if perturbation to an alternative state was 282 possible in this simple case (see Supplementary File 7). We found that only 0.079% of 3-283 species Lotka-Volterra systems exhibit the behaviour required by our two-state model, 284 suggesting that this model is unrealistic for small numbers of species (as we assume that 285 diversity is a continuous variable). However, for larger numbers of species, theoretical 286 ecology gives a strong justification for our assumptions. It has recently been shown that as 287 the number of species ݊ increases, the number of fixed points which are stable increases parameters that have multiple fixed points also increases: with ݊ ൌ 4 0 0 , this proportion is 290 >97% (36). This suggests that the overwhelming majority of mathematically possible 291 systems at relevant numbers of species exhibit multiple fixed points; the fraction of 292 biologically possible systems exhibiting this behaviour is likely even higher. Furthermore, 293 when resource competition is incorporated -a more realistic assumption in the case of 294 the human microbiome -all these fixed points become stable or marginally stable (36). 295 The gut microbiome is an ecosystem of hundreds of species in the presence of resource 296 competition. Goyal et al. recently showed that multiple resilient stable states can exist in 297 microbial communities if microbes utilize nutrients one at a time (22). We can therefore 298 state confidently that: the gut microbiome exists with multiple stable equilibria; its 299 community composition is history-dependent; and perturbations lead to transitions 300 between the multiple possible stable states. All of these assumptions justify the simplistic 301 coarse-grained model we describe here, which effectively takes these high-level emergent 302 properties of multi-species Lotka-Volterra models to build a substantially simpler model 303 based on a single, commonly-used metric: diversity. 304

305
Starting from a common conceptual picture of the microbiome as resting within a stability 306 landscape, we have developed a mathematical model of its response to perturbation by 307 antibiotics. Our framework, based on phylogenetic diversity, successfully captures the 308 dynamics of a previously published dataset for four common antibiotics (3), providing 309 quantitative support for these simplifying ecological assumptions. Using model selection, state transition in the oral microbiome with clindamycin, which was not detected by the 312 initial authors using a GLM repeated measures test. 313 While pairwise comparisons based on diversity can still identify differences in microbiome 314 state, they provide no information on microbiome dynamics. Our dynamical systems 315 approach therefore also gives additional mechanistic insight in this regard. Zaura et al. 316 observed that the lowest diversity in the gut microbiome was observed after a month rather 317 than immediately after treatment stopped (3). This cannot be due to a persistence of the 318 antibiotic effect, as all antibiotics used only have short half-lives of the order of hours 319 (37,38). Within our framework, this is because the full effects of the transient impulse take 320 time to be realised due to the overdamped nature of the system. We found a consistent 321 damping ratio for both ciprofloxacin and clindamycin, supporting this conclusion. 322 We have also demonstrated how our modelling framework could be used to compare 323 different hypotheses about the long-term effects of antibiotic perturbation by fitting different 324 models and using Bayesian model selection. Our modelling work provides an additional 325 line of evidence that while short-term restoration obeys a simple impulse response model, 326 the underlying long-term community state can be fundamentally altered by a brief course 327 of antibiotics, as suggested previously by others (7), raising concerns about the long-term 328 impact of antibiotic use on the gut microbiome. While this state transition may not 329 necessarily equate to any negative health impacts for the host (none of the participants 330 involved in the original study reported any gastrointestinal disturbance), in the gut 331 microbiome the transition to a new state with reduced diversity may increase the risk of 332 colonisation and overgrowth of pathogenic species. Interestingly, in the salivary 333 microbiome the transition appeared to be to a state with increased diversity, which is 334 associated with a greater risk of disease in the oral cavity (39). This observation was not noted by Zaura et al. -a significant difference was detected in diversity in one antibiotic 336 but was relegated to a supplementary figure and not discussed (3) -perhaps because it 337 appeared contradictory to their other conclusions. However, we believe it makes sense 338 within a stability landscape framework. Even if only marginal, when considered at a 339 population level these effects may mean that antibiotics have substantial negative health 340 consequences which could support reductions in the length of antibiotic courses, 341 independently of concerns about antibiotic resistance (40). Modelling the long-term impact 342 on the microbiome of different doses and courses could help to influence the use of 343 antibiotics in routine clinical care. Our sample size is small, so the precise posterior 344 estimates for parameters that we obtain should not be over-interpreted, but comparing 345 antibiotics using these parameter estimates represents another practical application. 346 Our framework lends itself naturally to comparing different dynamical models. We see our 347 two variant models as a starting point for a stability landscape approach, and would hope 348 that better models can be constructed. Hierarchical mixed effects models may offer an 349 improved fit, particularly if they take into account other covariates; however, we lacked the 350 necessary metadata on the participants from the original study (Table 1)  We would not expect the behaviour of the microbiome after longer or repeated courses of 363 antibiotics to be well-described by an impulse response model which assumes the course 364 is of negligible duration. Nevertheless, it would be possible to use the mathematical 365 framework given here to obtain an analytic form for the possible system response by 366 convolving any given perturbation function with the impulse response. It remains to be 367 seen whether this simple model would break down in such circumstances. 368 As we have demonstrated, while the individualized nature of the gut microbiome's 369 response to antibiotics can be highly variable, a general model still captures important 370 microbiome dynamics. We believe it would be a mistake to assume that our model is 'too 371 simple' to provide insight on a complex ecosystem. At this stage of our understanding, 372 creating a comprehensive inter-species model of the hundreds of members of the gut 373 microbiome appears intractable; it may also not be necessary for building simple models to 374 inform clinical treatment based on limited and sparse data. We believe there is a place for 375 both fine-grained models using pairwise interactions -particularly for systems of reduced 376 Fitting the model therefore requires fitting three parameters: Resulting in the following model (Model 1, Figure 1C): 415 Antibiotics may lead not just to displacement from equilibrium, but also state transitions to 417 new equilibria (2). To investigate this possibility, we also consider a model where the value 418 of equilibrium diversity asymptotically tends to a new value ‫ܣ‬ (Model 2, Figure 1C). As we 419 are aiming to minimise model complexity, we do this by adding a single parameter and a 420 term that asymptotically grows as time increases: 421  Files 3 and 4). 441 We found no association between sequencing depth and timepoint. 442

Phylogenetic diversity 443
There are many possible diversity metrics that could be used to compute the displacement 444 from equilibrium. Because of our assumption that phylogenetic diversity approximates assumptions'), we chose to use Faith's phylogenetic diversity (44)  We used a Bayesian framework to fit our basic model 1 (eq. 3) using Stan (48) and RStan 457 (49) to the gut and oral microbiome samples for the five separate groups: placebo, 458 ciprofloxacin, clindamycin, minocycline, and amoxicillin (i.e. n=2x5=10 fits). In brief, our 459 approach used 4 chains with a burn-in period of 1 000 iterations and 9 000 subsequent 460 iterations, verifying that all chains converged (ܴ = 1) and the effective sample size for each 461 parameter was sufficiently large (neff > 1 000). We additionally fitted model 2 with a 462 possible state transition (eq. 4) to all non-placebo groups (n=2x4=8 fits   Ciprofloxacin 9 9 Clindamycin 9 9 Minocycline 10 10 Amoxicillin 12 12