Social interaction, noise and antibiotic-mediated switches in the intestinal microbiota.

The intestinal microbiota plays important roles in digestion and resistance against entero-pathogens. As with other ecosystems, its species composition is resilient against small disturbances but strong perturbations such as antibiotics can affect the consortium dramatically. Antibiotic cessation does not necessarily restore pre-treatment conditions and disturbed microbiota are often susceptible to pathogen invasion. Here we propose a mathematical model to explain how antibiotic-mediated switches in the microbiota composition can result from simple social interactions between antibiotic-tolerant and antibiotic-sensitive bacterial groups. We build a two-species (e.g. two functional-groups) model and identify regions of domination by antibiotic-sensitive or antibiotic-tolerant bacteria, as well as a region of multistability where domination by either group is possible. Using a new framework that we derived from statistical physics, we calculate the duration of each microbiota composition state. This is shown to depend on the balance between random fluctuations in the bacterial densities and the strength of microbial interactions. The singular value decomposition of recent metagenomic data confirms our assumption of grouping microbes as antibiotic-tolerant or antibiotic-sensitive in response to a single antibiotic. Our methodology can be extended to multiple bacterial groups and thus it provides an ecological formalism to help interpret the present surge in microbiome data.

Introduction

Recent advances in metagenomics provide an unprecedented opportunity to investigate the intestinal microbiota and its role in human health and disease [1], [2]. The analysis of microflora composition has a great potential in diagnostics [3] and may lead to the rational design of new therapeutics that restore healthy microbial balance in patients [4]–[6]. Before the clinical translation of human microbiome biology is possible, we must seek to thoroughly understand the ecological processes governing microbiota composition dynamics and function.

The gastro-intestinal microbiota is a highly diverse bacterial community that performs an important digestive function and, at the same time, provides resistance against colonization by entero-pathogenic bacteria [7]–[9]. Commensal bacteria resist pathogens thanks to resources competition [1], [8], growth inhibition due to short-chain fatty acid production [10], killing with bacteriocins [11], [12] and immune responses stimulation [13], [14]. However, external challenges such as antibiotic therapies can harm the microbiota stability and make the host susceptible to pathogen colonization [15]–[20].

Despite its importance to human health, the basic ecology of the intestinal microbiota remains unclear. A recent large-scale cross-sectional study proposed that the intestinal microbiota variation in humans is stratified and fits into distinct enterotypes, which may determine how individuals respond to diet or drug intake [21]. Although there is an ongoing debate over the existence of discrete microbiome enterotypes [22], they could be explained by ecological theory as different states of an ecosystem [23]. Ecological theory can also explain how external factors, such as antibiotics, may lead to strong shifts in the microbial composition. A recent study that analyzed healthy adults undergoing consecutive administrations of the antibiotic ciprofloxacin, showed that the gut microbiota changes dramatically by losing key species and can take weeks to recover [24]. Longitudinal studies, such as this one, suggest that many microbial groups can have large and seemingly random density variations in the time-scale of weeks [25], [26]. The observation of multiple microbial states and the high temporal variability highlight the need for ecological frameworks that account for basic microbial interactions, as well as random fluctuations [27]–[29].

Here we propose a possible model to study how the intestinal microbiota responds to treatment with a single antibiotic. Our model expands on established ecological models and uses a minimal representation with two microbial groups [30] representing the antibiotic-sensitive and antibiotic-tolerant bacteria in the enteric consortium (Fig. 1). We propose a mechanism of direct interaction between the two bacterial groups that explains how domination by antibiotic-tolerants can persist even after antibiotic cessation. We then develop a new efficient framework that deals with non-conservative multi-stable field of forces and describes the role played by the noise in the process of recovery. We finally support our model by analyzing the temporal patterns of metagenomic data from the longitudinal study of Dethlefsen and Relman [24]. We show that the dynamics of microbiota can be qualitatively captured by our model and that the two-group representation is suitable for microbiota challenged by a single antibiotic. Our model can be extended to include multiple bacterial groups, which is necessary for a more general description of intestinal microbiota dynamics in response to multiple perturbations.

Figure 1

The two-group model of the intestinal microbiota with antibiotic-sensitive and antibiotic-tolerant bacteria.

Antibiotic sensitives can inhibit the growth of tolerants and both groups compete for the same growth substrate. Model parameters and represent the antibiotic sensitivity of sensitive and tolerant bacteria (where ), and represent their affinities to substrate and represents the inhibition of tolerants by sensitives.

Results

Mathematical model

We model the microbiota as a homogeneous system where we neglect spatial variation of antibiotic-sensitive () and antibiotic-tolerant () bacterial densities. Their evolution is determined by growth on a substrate and death due to natural mortality, antibiotic killing and social pressure. With these assumptions, we introduce, as a mathematical model, two coupled stochastic differential equations for the density of sensitives and tolerants ( and ) normalized with respect to the maximum achievable microbial density: In the physics literature these types of equations represent stochastic motion in a non-conservative force field . The first terms in correspond to the saturation growth terms representing the indirect competition for substrate and depend on , which is the ratio of the maximum specific growth rates between the two groups. If tolerants grow better than sensitives on the available substrate and the reverse is true for . They effectively describe a microbial system with a growth substrate modeled as a Monod kinetic [31] in the limit of quasi-steady state approximation for substrate and complete consumption from the microbes (see Methods for details). Both groups die with different susceptibility in response to the antibiotic treatment, which is assumed to be at steady-state. defines the ratio of the combined effect of antibiotic killing and natural mortality rates between the two groups (see Methods for details). While the system can be studied in its full generality for different choices of , we consider the case of because it represents the more relevant case where sensitives are more susceptible to die than tolerants in the presence of the antibiotic. A possible that mimics the antibiotic treatments is a pulse function. With this, we are able to reproduce realistic patterns of relative raise (fall) and fall (raise) of sensitives (tolerants) due to antibiotic treatment as we show in Fig. S4 in the Supporting Information Text (Text S1). Additionally, we introduce the social interaction term between the two groups, , to implement competitive growth inhibition [13], [32]. In particular, we are interested in the case where the sensitives can inhibit the growth of the tolerants (), which typically occurs through bacteriocin production [33]. Finally we add a stochastic term that models the effect of random fluctuations (noise), such as random microbial exposure, which we assume to be additive and Gaussian. The analysis can be generalized to other forms of noise such as multiplicative and coloured.

Antibiotic therapy produces multistability and hysteresis

We first analyzed the model in the limit of zero noise, . In this case, we were interested in studying the steady state solutions that correspond to the fixed-points of equations (1,2) and are obtained imposing . We found three qualitatively-distinct biologically meaningful states corresponding to sensitive domination, tolerant domination and sensitive-tolerant coexistence (see Text S1). We evaluated the stability of each fixed point (see Text S1) and identified three regions within the parameter space (Fig. 2A). In the first region the effect of antibiotics on sensitive bacteria is very low () and domination by sensitives is the only stable state (sensitives monostability). In the second region the effect of the antibiotic on sensitives is stronger than their inhibition over tolerants () and the only stable state is domination by tolerants (tolerants monostability). Finally, in the third region () both sensitive and tolerant dominations are possible and stable, while the third coexistence fixed point is unstable (bistability) (see Text S1). This simple analysis shows that multistability can occur in a gut microbiota challenged by an antibiotic where one group directly inhibits the other (i.e. through the term). Furthermore, it suggests that multistability is a general phenomenon since it requires only that antibiotic-sensitive and antibiotic-tolerant bacteria have similar affinities to nutrients. This is a realistic scenario because tolerants, such as vancomycin resistant Enterococcus [18], are often closely related to other commensal but antibiotic-sensitive strains and therefore should have similar affinity to nutrients. Finally, the solution of equations (1) and (2) reveals that hysteresis is present for values of and leading to multistability (Fig. 2B). Similarly to magnetic tapes, such as cassette or video tapes, which remain magnetized even after the external magnetic field is removed (i.e. stopping the recording), a transient dose of antibiotics can cause a microbiota switch that persists for long time even after antibiotic cessation.

Figure 2

Multistability and hysteresis in a simple model of the intestinal microbiota. A: phase diagram showing the three possible stability regions.

Antibiotic-sensitive bacteria dominate when and antibiotic-tolerant bacteria dominate when and therefore these are regions of monostability. There is a region of bistability between the two regions where domination by either sensitives or tolerants is possible. B: schematic display of the hysteresis phenomenon explaining cases where antibiotic treatment produces altered microbiota (i.e. tolerants domination) that persists long after antibiotic cessation. C–F: mean density values obtained simulating the Langevin dynamics for a maximum time after an instantaneous change of the parameter (C and D) and (E and F). These averages are obtained over noise realizations. C, D and E, F show the antibiotic-tolerants or antibiotic-sensitives densities, respectively, as a function of the social interaction parameter () with or the antibiotic killing () with .

Noise alters stability points

The previous analysis shows the existence of multistability in the absence of noise. However, the influx of microbes from the environment and/or the intra-population heterogeneity are expected in realistic scenarios and affect the bacterial density evolution in a non-deterministic fashion. This raises the question of how the noise alters the deterministic stable states and their stability criteria. We assume that the noise is a fraction of bacteria that can be added (or removed) at each time step, but on average has no effect since . This assumption is justified by the fact that a persistent net flux of non-culturable bacteria from the environment is unrealizable. We also assume that this random event at time is not correlated to any previous time , which corresponds to , where characterizes the noise amplitude and is the Dirac delta function. We calculated the stationary probability of the microbiota being at a given state by solving the stationary Fokker-Planck Equation (FPE) [34] corresponding to the Langevin equations (1,2):By numerically solving equation (3) as described in [35], for increasing , we find that for small values of the most probable states coincide with the deterministic stable states given by (Fig. 3A). However, by increasing the distribution spreads and the locations of the most probable states change and approach each other. As a consequence, the probability of an unstable coexistence, characterized by and , increases thus avoiding extinction. This intuitively justifies how recovery to a sensitive-dominated state within a finite time after antibiotic cessation becomes possible with the addition of the noise. Without noise, the complete extinction of sensitive bacteria would have prevented any possible recolonization of the intestine. Beyond a critical noise level () bistability is entirely lost and the probability distribution becomes single-peaked with both bacterial groups coexisting. The microbiota composition at the coexistence state can be numerically determined from the solution of , as shown in Fig. 3B and Video S1. Further investigations based on analytical expansion of the Langevin equations (see Methods) show that for small random fluctuations, , the first noise-induced corrections to the deterministic density are linearly dependent on with a proportionality coefficient determined by the nature of the interactions (insets in Fig. 3B). These linear correction terms can be obtained as a function of the model parameters and, after substituting a particular set of values in the bistable region ( and ), they are for sensitives and for tolerants. These numbers are different from those reported in the insets of Fig. 3B. However the discrepancy is due to the propagation of the boundary conditions when numerically solving the solution of the FPE using finite elements (see Text S1).

Figure 3

Most probable microbiota states change from bistable scenario to mono-stable coexistence with increasing noise.

A: the bacterial density joint probability distribution determined by solving the Fokker-Planck equation (3) for four different values of the environmental noise. B: the bacterial densities at the peaks of as a function of the noise parameter . Red symbols are data from the numeric solution of the Fokker-Planck equation and the black solid lines are the exponential fit. Parameters used: and . The insets detail the linear regime.

This has important biological implications since it suggests that extinction is prevented and, more importantly, that a minority of environmental microbes can settle in the gut at a rate that depends on the strength of their social interaction with the established microbiota.

The introduction of random perturbation affects the stability criteria of the stable states. In particular, we observe that the bistability region decreases when the noise amplitude increases (Fig. 2C–F). At the limit, when the bistability is entirely lost and the only stable state is the one where both groups coexist. This concept was previously hypothesized but not explicitly demonstrated in a model of microbial symbionts in corals [30].

Noise affects the recovery time

Our model predicts that in absence of stochastic fluctuations the recovery time is larger than any observational time-scale so that it is impossible to revert to the conditions preceding antibiotic perturbation (see Fig. S4 in Text S1). In reality, data show that this time can be finite and depends on the microbiota composition and the degree of isolation of the individuals [18], [24], [36]. Thus, we aim to quantitatively characterize how the relative contribution of social interaction and noise level affects the computation of the mean residence time.

In order to determine the relative time spent in each domination state, we compute the probability of residence in each stable state using master equations [34]. This method is more efficient than simulating the system time evolution by direct integration of the Langevin equations because it boils down to solving a deterministic second-order differential equation. Furthermore, this approach scales up well when the number of microbial groups increases, in contrast to the numerical solution of the FPE which can become prohibitive when . In our model, the master equations for the probability of residing in the tolerant or sensitive domination state are: where is the transition rate from state to , which can be obtained in terms of the sum over all the state space trajectories connecting to .

By solving this system of equations at steady-state, we obtain the residence probabilities and . After computing the transition rate as a function of the parameters, as reported in the Methods, we determine , which is our theoretical prediction for the mean relative residence time spent in the tolerant domination state (Fig. 4). The theoretical predictions are in good agreement with those obtained by simulating the dynamics multiple times and averaging over different realizations of the noise. A first consequence from this analysis is that the time needed to naturally revert from the altered state depends exponentially on the noise amplitude (). As such, we predict that for the case of an isolated system () the switching time is exponentially larger than any other microscopic scale and the return to a previous unperturbed state is very unlikely. On the contrary, as the level of random exposure is increased, the time to recover to the pre-treated configuration decreases (see Fig. S4 in Text S1). Additionally, this method can be considered as a way to indirectly determine the strength of the ecological interactions between microbes which can be achieved by measuring the amount of time that the microbial population spends in one of the particular microbiota states. Therefore, it can potentially be applied to validate proposed models of ecological interactions by comparing residence times measured experimentally with theoretical predictions.

Figure 4

Microbiota resident time in antibiotic-tolerant domination as a function of the: A–B) antibiotic action () and C–D) social interaction () parameters.

Blue circles show the theoretical predictions obtained by determining the probability of the most probable path. Red circles are obtained by simulating the Langevin dynamics over iterations and averaged for noise realizations. Higher order-corrections can be included to increase the theoretical estimation accuracy.

Analysis of metagenomic data reveals antibiotic-tolerant and antibiotic-sensitive bacteria

We now focus on the dynamics of bacteria detected in the human intestine and test the suitability of our two-group representation by re-analyzing the time behaviour in the recently published metagenomic data of Dethlefsen and Relman [24]. The data consisted of three individuals monitored over a 10 month period who were subjected to two courses of the antibiotic ciprofloxacin. Since the data are noisy and complex, and the individual subjects' responses to the antibiotic are distinct [24], identifying a time behaviour by manual screening is not a trivial task. We do it by using singular value decomposition (SVD) to classify each subject phylotype-by-sample data matrix into its principal components. Because of inter-individual variability we obtain, for each subject, the right and left eigenvectors associated to each eigenvalue. By ranking the phylotypes based on their correlation with the first two components we recover characteristic temporal patterns for each volunteer [37], [38].

In all three subjects, we observe that, in spite of the individualized antibiotic effect, the two dominant eigenvalues or principal components together capture about 70% of the variance observed in the data (Fig. 5A–C). Invariably, the first component shows a decrease in correspondence to antibiotic treatment and reflects the behaviour of antibiotic-sensitive bacteria (green line in Fig. 5D–F). Conversely, the second component increases with the antibiotic treatment and represents antibiotic-tolerants (red line in Fig. 5D–F). The observation that each subject's microbiota can be decomposed into two groups of bacteria with opposite responses to antibiotics supports the validity of the two-group approach used in our model. Classification of each individual's phylotypes as sensitive or tolerant can be obtained by determining their correlation with the two principal components (see Text S1) (information in the right-eigenarrays matrix from SVD). Bacteria correlated with component 1 are usually highly abundant before antibiotic treatment and drop strongly during treatment, often below detection. Vice-versa, bacteria correlated with component 2 are typically in low abundance before the antibiotic and increase with antibiotic administration (Fig. 5G–I). Interestingly, despite significant inter-individual differences in recovery time (Fig. 5G–I) and individualized response of each subject, the data show that in each individual the majority of bacteria are antibiotic-sensitive and only a small but significant fraction are tolerant to ciprofloxacin (see Text S1). The recognition of these time-patterns could be considered as a possible tool to indirectly determine the susceptibility of non-culturable commensal bacteria to FDA-approved antimicrobial compounds. However, the presence of strains in the same phylotypes that display both behaviors in response to the drug may constitute a significant challenge for the success of this method.

Figure 5

Analysis of microbiota response to the antibiotic ciprofloxacin from three subjects [24] using singular value decomposition identifies antibiotic-sensitive and antibiotic-tolerant bacteria.

A–C: fraction of variance explained by the five most dominant components. D–F: plot of each sample component 1 (green) and 2 (red) coordinates versus sample time. G–I: sorting of the phylotypes log2-transformed abundance matrix based on the correlation within the two principal component. Above (below) the green dashed lines, we display the time series of the top 20 phylotypes strongly correlated (anti-correlated) with component 1 and anti-correlated (correlated) with 2 and dropping (increasing) during treatment, which we identify as sensitves (tolerants). Subject 3 (C,F,I) displays absence of sensitive bacteria for a prolonged period of about 50 days after the first antibiotic treatment. This confirms the fact that microbiota response to antibiotic can differ from subject to subject. Additionally, it also supports our model prediction of remaining locked in a tolerant-dominated state after antibiotic treatment cessation.

The time evolution of the phylotypes (Fig. 5G–I) qualitatively agrees with our theoretical prediction that after the antibiotic administration the system moves fast, meaning in a time smaller than any other observable time-scale, into a new stable state with less sensitives and more tolerants. Further, the data also suggest that the return to sensitive domination happens after a recovery-time scale that depends on the microbial composition.

Discussion

We present a model of inter-bacterial interactions that explains the effect of antibiotics and the counter-intuitive observation that an antibiotic-induced shift in microbiota composition can persist even after antibiotic cessation. Our analysis predicts a crucial dependence of the recovery time on the level of noise, as suggested by experiments with mice where the recovery depends on the exposure to mice with untreated microbiota [18]. The simple model here introduced is inspired by classical ecological modeling such as competitive Lotka-Volterra models [39], [40], but relies on mechanistic rather than phenomenological assumptions, such as the logistic growth. Although more sophisticated multi-species models include explicit spatial structure to describe microbial consortia [33], [41]–[43], our model is a first attempt to quantitatively analyze the interplay between microbial social interactions () and stochastic fluctuations () in the gut microbiota. We find that these two mechanisms are the key ingredients to reproduce the main features of the dynamics in response to antibiotic (sudden shifts and recovery). Our model can be easily generalized to include spatial variability and more complicated types of noise. Therefore we provide a theoretical framework to quantify microbiota resilience against disturbances, which is an importance feature in all ecosystems [44]. By introducing a new stochastic formulation, we were able to characterize composition switches within the context of state transition theory [45], [46], an important development over similar ecological models of microbial populations [30]. We present a new method to calculate the rate of switching between states that identifies the most likely trajectory between two stable states and their relative residence time, which can be subjected to experimental validation. Finally, we apply SVD to previously published metagenomic data [24], which allows us to classify the bacteria of each subject in two groups according to their temporal response to a single antibiotic. The SVD method has been used before to find patterns in temporal high-throughput data, including transcription microarrays [37] and metabolomics [47]. Although our approach seems to capture well the main temporal microbiota patterns, we should note that the use of the Euclidean distance as a metric for microbiome analysis presents limitations and recent studies have proposed alternative choices [48]–[50]. We also opt for an indirect gradient analysis method [51] because we are interested in emergent patterns from the data regardless of the measurements of the external environmental variable (i.e. presence or absence of the antibiotic) [50].

We propose a mechanism of interaction between two bacterial groups to explain the lack of recovery observed in the experiments that can be validated in the near future. Although training the model with the available data sets would be of great interest, this will not be useful in practice because we need more statistical power to be predictive. However, we anticipate that a properly validated mathematical model of the intestinal microbiota will be a valuable tool to assist in the rational design of antibiotic therapies. For example, we predict that the rate of antibiotic dosage will play a crucial role. In order to let the microbiota recover from antibiotic treatment, it is better to gradually decrease antibiotic dosage at the first sign of average microbiota composition change, which has to be larger than the threshold community change represented by the day-to-day variability [26], rather than waiting for tolerant-domination and then stopping antibiotic treatment.

We show here the application of our theory to a two-bacterial group scenario because we are interested in the microbiota response when challenged with a single antibiotic. However, in more realistic conditions the microbiota is subjected to different types of perturbations, which may drive it towards more alternative stable states. Our theory of the microbial-states switches characterization can be naturally extended to more than two states and consists of the solution of the linear system of equations , where is the array of probability of residing in each stable state and is the matrix of transition rates among the states.

The ongoing efforts to characterize the microbial consortia of the human microbiome can yield tremendous benefits to human health [52]–[55]. Within the next few years, we are certain to witness important breakthroughs, including an increase in the number of microbiomes sequenced as well as in sequencing depth. Yet, without the proper ecological framework these complex ecosystems will remain poorly understood. Our study shows that, as in other complex microbial ecosystems, ecological models can be valuable tools to interpret the dynamics in the intestinal microbiota.

Methods

Full model and simplification

The model introduced in equations 1 and 2 is derived from the more detailed model described below. We model the bacterial competition in a well-mixed system in the presence of antibiotic treatment by means of the following stochastic differential equations: where we account for two bacterial groups; the intestinal resident sensitive flora and an antibiotic tolerant one . Additionally, we also consider the substrate and the antibiotic densities. The antibiotic time evolution is simply a balance between inflow and outflow (i.e. no decay due to microbial degradation) where is the system's dilution rate, which sets the characteristic microscopic time-scale, and is the constant density of the incoming antibiotic, which can be time dependent. Similarly the substrate concentration, , results from a mass balance from influx and microbial consumption. As for the antibiotic, is the constant density of the incoming nutrient (i.e. the concentration of resources coming from the small-intestine). The second and third terms in the right-hand side of the second equation in (5) describe the amount of substrate consumed by bacterial growth assuming Monod kinetics where () is the maximum growth rate for sensitives (tolerants), is the half-saturation constant for growth, which parametrizes the bacterial affinity to the nutrient, and () is the yield for growth for sensitives (tolerants). The last two equations describe how sensitives and tolerants grow on the substrate available and are diluted with the factor . We mimic the effect of the antibiotic on the sensitives adding a term proportional to the sensitive density where the constant of proportionality is the antibiotic-killing rate. We also introduce a direct inhibition term , which mimics the inhibition of sensitive bacteria on the tolerants (social interaction). Finally the Gaussian random variables , are the additive random patterns of exposure and represent the random microbial inflows (outflows) from (to) the external environment.

It is convenient to scale the variables and set the dilution rate to unity (). Therefore, all the rates have to be compared with respect to the system characteristic dilution rate. Introducing , , , , , , , , , and and dropping the tilde symbols, we obtain the following dimensionless model: If we assume that the antibiotic is a fast variable compared to the microbial densities () (i.e. the time-scale at which the antibiotic reaches stationary state is smaller than that of the bacteria), we can solve for and obtain . If we also assume that the incoming substrate is all consumed in microbial growth, therefore maintaining the population in a stationary state with respect to the available resources, and that, similarly to the antibiotic, the resources equilibrate much faster than the bacterial densities (quasi-steady state assumption, ), we obtain that:If we now define a new parameter describing the relative ratio of the combination of antibiotic killing and natural mortality (i.e. wash-out) between sensitives and tolerants, the model reduces to the two variables model in reported in equations (1–2).

Effective potential and location of long-term states

The introduction of random noise has the important consequence of changing the composition of the stable states (Fig. 3A). In order to characterize this phenomenon, we expand the solution of the Langevin equations (1–2) around one of the stable states obtaining the following set of equations for the variable :where to simplify the notation we drop the explicit time-dependence. We can easily recognize the first derivative of the force on the right-hand side as the Jacobian matrix computed in one of the minima . This equation can be solved order by order by defining the expansion and writing the equations for each order as:

Assuming that the initial condition at time zero is , which can always be neglected for long-term behaviour, the solution of equation (9) isThis means that the average location of the minima at zero order is not modified by the noise since . By computing the solution of the equation (10) we similarly find that:The long-time average value of the first order correction now reads:The time integral can be easily computed assuming that the eigenvalues of are negative, or at least their real part is, as it should be for stable fixed points; therefore we obtain that:Thus, we find that the effect of random fluctuations is to correct the value of the stable points as if an external field, proportional to strength of the fluctuations, was present. This field is equal to the mean square displacement at large time opportunely weighted by the inverse of the curvature of the bare potential around the stable points, . The correlation can be now computed using equation (11) and reads:Since the previous equation simplifies towhich results in .

Theoretical estimate of the mean residence time

The mean residence time in each state is proportional to the residence probability defined in equation (4). To obtain it, we need to compute the transition rate as a function of the model parameters as:where and are the initial and final time and is the functional integral over the trajectory . Each time trajectory , solution of equations (1–2), has an associated weight , defined as:By discretizing the time so that with and the microscopic time step, we obtain that the Langevin equations can be written using the Ito prescription [56] as:where we use the short notation and the initial value is . The time discretization allows us to interpret the functional integral in equation (18) as:Since the noise is Gaussian and white, its distribution now reads:This can be justified using the property of the delta-function and its discrete time version so that follows and .

Using the properties of the delta function, and integrating out all s, the continuous limit expression of equation (21) iswhere has an intuitive interpretation in thermodynamics and it is related to the entropy production rate [57]. By using stationary-phase approximation, it turns out that in the computation of the rate defined in (17) only one path matters, , which is the most probable path. Higher order factors are proportional to the term [45], [46], and therefore simplify with the denominator in equation (21). This comes from the fact that several almost optimal paths can be constructed starting from . In the optimal path, the system stays in a stable state for a very long time, then it rapidly switches to the other stable state where it persists until . By shifting the switching time one obtains sub-optimal paths that, at the leading order in , give the same contribution of the optimal one and their number is directly proportional to . This leads to

The functional Gaussian integral can be computed [45], [46] and only provides a sub-leading correction to the saddle-point contribution resulting in the transition rate formula , which is reported in the Results section.

We now need to determine the optimal path and its associated action . This path is defined as the one where the functional derivative of is set to zero such that the initial and final states are fixed. This produces a set of second-order differential equationswhich can be solved imposing the initial conditions on and .

It is easy to verify that the downhill solution is and it is associated with null action. Meanwhile, the ascending trajectory, which is the one leading to a non-zero action and hence gives the transition rate value, is not given by , as it would be for conservative field of forces. This means that in presence of a dissipative term the reverse optimal path from the minimum to the maximum is different with respect to the one connecting the maximum from the minimum of the landscape.

As the last point, we want to show that the action associated to the optimal path can be further simplified by noticing thatWe can easily prove this condition by showing that the time derivative vanishes when equation (24) is satisfied and remembering that the optimal path connects two stable states where and . This property allows us to rewrite the action as:

We solved numerically the equation (24) using a trial-and-error approach. We varied the first-derivative at initial time in order to arrive as close as possible to the final point within some numerical precision. In principle the ideal trajectory connecting two stable points should be computed in the limit of but this trajectory will take infinite time. We report three examples of most probable paths connecting the points to and reverse for a chosen set of in Fig. S6 of the Text S1.

Singular value decomposition

We first rarefy the raw phylotypes counts matrix as in [24]. We then normalize the logarithm of the counts according to the following procedure: 1) we add one to all the phylotypes counts to take into account also for the non-detected phylotypes in each sample, 2) we log-transform the data and 3) we normalize the resulting matrix with respect to the samples averages. In formulae, the count associated to phylotype in sample for each subject iswhere is the average value of the counts in each sample and is the total number of phylotypes. Among all possible normalization schemes, we decide to subtract the column averages because we aim at identifying patterns within samples based on their correlation in bacterial composition. Indeed, the covariance matrix of the samples is proportional to , where is the transpose matrix. SVD on the matrix is thus equivalent to the principal component analysis (PCA) performed on the samples covariance matrix.

Text S1 reports additional calculations, figures and details on: 1) model and relative stability analysis, 2) effect of random fluctuations and noise-induced dynamics and 3) Singular Value Decomposition.

(PDF)

Click here for additional data file.

Video S1 shows the stationary probability distributions as a function of the sensitive and tolerant densities for increasing noise value , which ranges from to . For visualization purposes, the noise value associated to each movie frame is displayed as an increasing bar in the top panel.

(MOV)

Click here for additional data file.

Video S2 shows the time evolution of the two principal components for the three subjects from [24]. Empty circles represent untreated samples, asterisks represent samples during treatment and filled circles represent represent samples during treatment 2.

(MP4)

Click here for additional data file.