Navigation and control of outcomes in a generalized Lotka-Volterra model of the microbiome

In this chapter the generalized Lotka-Volterra (gLV) equations are treated as a mathematical proxy for the microbial dynamics of the gut microbiome. The gLV equations describe the dynamics of a microbial population $y_{i}$ in a system of $N$ total interacting populations as

$$\frac{dy_{i}}{dt}=y_{i}\left(\rho_{i}+\sum_{j=1}^{N}K_{ij}y_{j}\right),\quad i\in 1,\ \ldots,\ N,$$

(1)

where the growth rate of species $i$ is given by $\rho_{i}$, and the pairwise interaction $K_{ij}$ encodes the ecological effects of species $j$ on species $i$.

A gLV system with $N$ populations can exhibit up to $2^{N}$ steady states, where each steady state is specific to a distinct presence/absence combination of the $N$ species. Motivated by clinical bacteriotherapies like fecal microbiota transplantation (FMT) that seek to drive a diseased microbiome towards a healthy composition, we consider how to control the steady-state outcomes of gLV systems with direct interventions that instantaneously supplement or deplete the microbial state of a system. Towards this end, as in previous work stein2013ecological ; jones2018silico the gLV model is extended to include the instantaneous addition of a foreign microbial transplant v at time $t^{*}$, as well as the administration of an antibiotic treatment $u(t)$, so that

$$\frac{dy_{i}}{dt}=y_{i}\left(\rho_{i}+\sum_{j=1}^{N}K_{ij}y_{j}\right)+v_{i}\delta(t-t^{*})+u(t)\varepsilon_{i}y_{i},$$

(2)

where $v_{i}$ is the $i$th component of the foreign transplant v, $\delta(\cdot)$ is the Dirac delta function, and $\varepsilon_{i}$ is the antibiotic susceptibility of population $i$.

Three types of bacteriotherapy-inspired direct interventions are examined, corresponding to three types of transplant v: the addition of a single species (interpreted as a probiotic), the addition of a stable external microbial community (interpreted as fecal microbial transplant), and the removal of a particular species (interpreted as a phage therapy). More explicitly, we consider the direct control problem in which an initial condition is in the basin of attraction of some attractor $\mathbf{y_{a}}$, and the goal is to identify the smallest intervention $\mathbf{v}$ that is able to drive the system into the basin of attraction of some other target state $\mathbf{y_{b}}$.

We choose the cost function $\varphi$ to be the size of the transplant,

$$\varphi=\lVert\textbf{v}\rVert_{1},$$

(3)

so that the optimal control protocol minimizes the required transplant size. Since the transplant inputs are instantaneous (i.e. Dirac-delta type), this is a sampled-data control problem BourdinTrelat2016 .

We use a gLV model of C. difficile infection (CDI) as a case study for how to transition between basins of attraction. In this model, the growth rates $\bm{\rho}$ and interaction parameters $\mathbf{K}$ were fit by Stein et al. to microbial abundance time-series data from a mouse experiment performed by Buffie et al. stein2013ecological ; BuffiePamer2012 . To reduce the number of parameters of the model, Stein et al. assumed that bacteria within a given genus behave similarly, and consolidated the species-level data into genus-level data. In Fig. 1 the interaction parameters are displayed as a food web where the circles are microbial populations, and the edges describe the effect of one population on another, where positive interactions ($K_{ij}>0$) are black and inhibitory interactions ($K_{ji}<0$) are red. The interactions between populations have no clear hierarchy, indicating that the fitted model describes microbes on the same trophic level competing for shared resources. The fitted antibiotic susceptibilities $\varepsilon_{i}$ are mostly negative, indicating that antibiotics tend to deplete the growth rates microbial populations.

The CDI model produces microbial trajectories that allow for the simulated application of medical interventions. Previous work explored the effects of antibiotics on microbial dynamics stein2013ecological ; jones2018silico in this CDI model, and found that the CDI model exhibits the clinically- and experimentally-observed outcome that antibiotic-treated microbiomes were vulnerable to CDI.

To demonstrate antibiotic-induced CDI in the gLV model, Fig. 2 shows simulated microbial trajectories that result from applying four separate intervention scenarios to an initial condition measured by Buffie et al. For clarity, these figures plot the total microbe count on a log scale (in which the total microbe count is the sum of all of the microbes in each microbial population), and at each time each microbial population is linearly colored according to its proportion at that time. First, in Fig. 2a the system is not perturbed and the system evolves to a “healthy” steady state, in the sense that CD is unable to invade this steady state. In Fig. 2b the system is exposed to a unit dose of antibiotics (i.e., $u(t)=1$ for $0\leq t\leq 1$, and $u(t)=0$ otherwise) which drives the system to an “antibiotic-depleted” state, in the sense that CD is able to invade it. Bearing this out, in Fig. 2c the system is exposed to a unit dose of antibiotics and then inoculated with CD on day 10 (i.e. the transplant v is a unit vector of CD applied at time $t^{*}=10$), and this system evolves towards an “CD infected” steady state in which CD is present. Finally, in Fig. 2d the system is exposed to a unit dose of antibiotics, inoculated with CD on day 1, and then also supplemented with an external foreign transplant (i.e. the foreign population v is composed of an experimentally-measured initial condition introduced at a time $t^{*}=1$); due to the direct intervention, this system now evolves towards the healthy steady state.

In each panel of Fig. 2 the colored mice represent attained steady-state microbiome compositions: green represents healthy, yellow represents antibiotic-depleted, and red represents CD-infected. These three steady states of the CDI model resemble the compositions of the experimentally-observed mouse microbiome compositions, including their susceptibility or resilience to CD exposure BuffiePamer2012 ; stein2013ecological . In addition to these three steady states, two other steady states can be reached by initializing the system at one of the nine experimentally-measured initial conditions and then applying some type of intervention. We call this set of five steady states the “reachable” steady states of the CDI model, and focus on the ecological dynamics near them. The microbial compositions of these five steady states are displayed in Fig. 3. In this figure, the previously-mentioned healthy state is labeled steady state C, the antibiotic-depleted state is labeled steady state E, and the CD-infected state is labeled steady state D.

This experimentally-characterized CDI model provides a clinically-motivated case study in which different steady states can be associated with biologically meaningful microbiome compositions. It is pertinent to be able to efficiently switch basins of attraction in order to attain a “healthy” state, and in the remainder of this chapter we investigate this goal in detail.

The ecological dynamics of high-dimensional gLV systems are straightforward to simulate, but difficult to investigate analytically. Often the most fundamental features of a model are captured not by the particular high-dimensional vector that describes the state of a system, but rather by an abstract biological outcome associated with the state of a system. For example, in the case of the CDI model, the most important question is whether the system will tend towards a healthy (steady state C) or an antibiotic-depleted state (steady state E), and the precise composition of those steady states is not as important. Thus, in some sense there is a low-dimensional description of the transition between a pair of steady states, even while the actual ecological dynamics flow in a high-dimensional state space.

In earlier work we exploited this abstract outcome-oriented perspective to design steady-state reduction (SSR) JonesCarlson2019 . This method compresses a bistable region of a high-dimensional gLV model into a reduced two-state gLV model in which the two basis populations correspond to the bistable states of the original model. As depicted in Fig. 4, this reduced two-dimensional (2D) gLV model is defined on the 2D subspace spanned by a pair of steady states of the original model, and the subspace itself is embedded within the high-dimensional ecological phase space of the original gLV model. The parameters of the reduced model are weighted combinations of the parameters of the original model, with weights that are related to the composition of the two high-dimensional steady states. Within this subspace, these reduced dynamics constitute the best possible 2D gLV approximation of the high-dimensional gLV dynamics. Though gLV systems in general are capable of displaying periodic or chaotic behaviors, here our attention is restricted to trajectories that approach a fixed point, and in particular to regions of phase space that are near the border of the basins of attraction of a pair of steady states (i.e., near the separatrix).

More explicitly, to determine the SSR-reduced 2D gLV system associated with a bistable high-dimensional gLV model, first define variables $x_{a}$ and $x_{b}$ in the direction of unit vectors $\hat{x}_{a}$ and $\hat{x}_{b}$ that parallel the two steady states according to $\hat{x}_{a}\equiv\vec{y}_{a}/\lVert\vec{y}_{a}\rVert_{2}$, and $\hat{x}_{b}\equiv\vec{y}_{b}/\lVert\vec{y}_{b}\rVert_{2}$, where $\lVert\cdot\rVert_{2}$ is the $2$-norm. The 2D gLV dynamics on the subspace spanned by $\hat{x}_{a}$ and $\hat{x}_{b}$ are given by

$\displaystyle\begin{split}\frac{\text{d}x_{a}}{\text{d}t}&=x_{a}(\mu_{a}+M_{aa}x_{a}+M_{ab}x_{b}),\quad\text{and}\\ \frac{\text{d}x_{b}}{\text{d}t}&=x_{b}(\mu_{b}+M_{ba}x_{a}+M_{bb}x_{b}).\end{split}$

(4)

Here, the in-plane dynamics on this subspace in vector form are written $\frac{\text{d}\vec{x}}{\text{d}t}=\frac{\text{d}x_{a}}{\text{d}t}\hat{x}_{a}+\frac{\text{d}x_{b}}{\text{d}t}\hat{x}_{b}$.

The 2D parameters generated by SSR are chosen to minimize the deviation between the in-plane and high-dimensional gLV dynamics $\epsilon=\lVert\frac{\text{d}\vec{y}}{\text{d}t}-\frac{\text{d}\vec{x}}{\text{d}t}\rVert_{2}$ for any point on the subspace spanned by $\hat{x}_{a}$ and $\hat{x}_{b}$. The values of the SSR parameters that minimize $\epsilon$ are derived in JonesCarlson2019 . When the two steady states $\vec{y}_{a}$ and $\vec{y}_{b}$ are orthogonal, the 2D parameters are given by

	$\displaystyle\mu_{\gamma}$	$\displaystyle=\frac{\vec{\rho}\cdot\vec{y}_{\gamma}^{\circ 2}}{\left\lVert\vec{y}_{\gamma}\right\rVert^{2}_{2}},\text{ and}$
	$\displaystyle M_{\gamma\delta}$	$\displaystyle=\frac{(\vec{y}_{\gamma}^{\circ 2})^{T}K\vec{y}_{\delta}}{\left\lVert\vec{y}_{\gamma}\right\rVert^{2}_{2}\left\lVert\vec{y}_{\delta}\right\rVert_{2}},$			(5)

where $\vec{y}^{\circ 2}\equiv\text{diag}(\vec{y})\vec{y}$ is the element-wise square of $\vec{y}$. When $\vec{y}_{a}$ and $\vec{y}_{b}$ are not orthogonal the interspecies interaction terms become slightly more complicated and are given by

	$\displaystyle M_{ab}$	$\displaystyle=\frac{\sum_{i,j=1}^{N}K_{ij}(y_{ai}y_{bj}+y_{bi}y_{aj})(y_{ai}-y_{bi}\sum_{k=1}^{N}y_{ak}y_{bk})}{1-(\sum_{i=1}^{N}y_{ai}y_{bi})^{2}},\text{ and}$
	$\displaystyle M_{ba}$	$\displaystyle=\frac{\sum_{i,j=1}^{N}K_{ij}(y_{bi}y_{aj}+y_{ai}y_{bj})(y_{bi}-y_{ai}\sum_{k=1}^{N}y_{bk}y_{ak})}{1-(\sum_{i=1}^{N}y_{bi}y_{ai})^{2}},$			(6)

where $\gamma,\delta\in a,b$, and $y_{ai}$ and $y_{bi}$ are the $i$th components of the of the unit vectors $\hat{y}_{a}\equiv\vec{y}_{a}/\lVert\vec{y}_{a}\rVert_{2}$ and $\hat{y}_{b}\equiv\vec{y}_{b}/\lVert\vec{y}_{b}\rVert_{2}$, respectively. Under this construction, the high-dimensional steady states $\vec{y}_{a}$ and $\vec{y}_{b}$ have in-plane steady state counterparts at $(\lVert\vec{y}_{a}\rVert_{2},\ 0)$ and $(0,\ \lVert\vec{y}_{b}\rVert_{2})$, respectively.

Crucially, this reduced 2D gLV system is analytically tractable: the separatrix can be written analytically JonesCarlson2019 , and bifurcation analyses readily inform how interaction parameters alter the basins of attraction of a system WangCarlson2020 . SSR associates complex high-dimensional dynamics with an intuitive low-dimensional system. In the context of the CDI model, this SSR-based understanding informs the transition between e.g. the healthy “green” and antibiotic-depleted “yellow” microbial states, and eschews the high-dimensional details. Later, we will decompose multistable gLV systems into a web of bistable subsystems, then use this web to describe the geometry of the basins of attraction of steady states of the model as an attractor network. Additionally, SSR will be leveraged to analytically compute the location of the separatrix in each subsystem, which is much more efficient than computing the locations of the separatrices numerically.

It is difficult to discern the influence of transplant size and timing on transplant efficacy in high-dimensional gLV systems, since only numerical methods are available to probe this relationship. However, SSR provides a link between high-dimensional and reduced 2D systems, allowing the mechanisms that underlie direct interventions in high-dimensional systems to be understood in terms of their low-dimensional counterparts.

Fig. 5 demonstrates how variability in transplant size and timing affects steady-state outcomes in a clinically-motivated scenario of CDI in a 2D gLV system. Here, an initial condition is depleted by antibiotics (RX) to the extent that it enters a diseased basin of attraction. Immediately, the reduced system makes clear that the location of the separatrix is crucial in determining the outcome of a microbial state: successful direct interventions must cross the separatrix. The transplant size required for success depends on the composition of the transplant (in this figure transplants are composed entirely of $\tilde{x}_{b}$). Furthermore, the required transplant size is variable, and depends on the ecological dynamics of the system— in Fig. 5 a transplant of the same size and composition is administered at two different times, but the later transplant is unsuccessful.

SSR allows for the relationships between transplant size, composition, and timing to be examined analytically in terms of an intuitive low-dimensional system. Next, having characterized the operation of direct transplants in a two-dimensional system, we examine how direct transplants may be used to transition between steady states in the full context of a multistable gLV system.

Having demonstrated that transplants are capable of shifting an ecological state into a different basin of attraction, we now investigate how transplant composition influences transplant efficacy. Three intervention types are considered: the introduction of a single-species “probiotic,” a stable community via “fecal microbiota tranplantation” (FMT), or the elimination of a single species via “phage therapy.” These interventions act as in silico proxies for medical therapies. The success rates of these interventions are plotted as a function of the intervention magnitudes in Fig. 6. As an example, the transplant administered in Fig. 2d is considered a “successful” intervention, since it was able to alter a trajectory tending towards steady state D (the CD-infected steady state) and drive it towards steady state C (the healthy steady state).

Fecal microbiota transplantation is implemented in the model by setting the transplant v proportional to a steady state of the gLV system. In particular, transplants v are set proportional to one of the five reachable steady states depicted in Fig. 3, or they are set proportional to one of the nine experimentally-measured initial conditions from the mouse experiment performed by Buffie et al. BuffiePamer2012 . Probiotics are realized by setting the transplant v proportional to a single microbial species. Lastly, phage therapies are described by making the transplant v negative, and setting it proportional to a single microbial species (so that it depletes a particular species).

In Fig. 6 these interventions are applied to initial conditions located at the reachable steady states of the CDI model depicted in Fig. 3. Interventions are considered successful when they shift the basin of attraction an initial condition is in. For the single species “probiotic,” transplants are solely composed of one of the 11 bacterial species of the model. The steady state “FMT” populations are composed of one of the five reachable steady states. There are 11 “phage” interventions that each deplete a single bacterial species. Lastly, the nine “experimental IC” foreign populations are proportional to the experimental initial conditions measured by Buffie et al. BuffiePamer2012 . In Fig. 6 the success rates of each intervention are plotted according to their magnitude (i.e. the one-norm $\lVert\textbf{v}\rVert_{1}$ of each foreign population). For scale, the five reachable steady states of the CDI model vary in size between $3\times 10^{11}$ and $24\times 10^{11}$ microbes, which informs the range of direct intervention sizes considered here.

As shown in Fig. 6, the success rates of multi-species interventions (steady states and experimental ICs) are significantly higher than the single-species interventions (single-species and phage therapies). In particular, phage therapies are completely ineffective at altering the basin of attraction that an initial condition is in. For each type of intervention, the success rates of each candidate transplant within each intervention type are plotted in dashed lines. The bulk success rates of each intervention are plotted in bold, and are computed by averaging the success rates of the individual introduced foreign populations within each intervention type.

The ecological dynamics of gLV systems are dictated by complex feedbacks between populations. Conceptually, in the absence of periodic or chaotic dynamics these ecological dynamics can be interpreted as flowing downhill a dynamical landscape (e.g. a Lyapunov function) towards a point attractor. A visualization of a dynamical landscape for a two-dimensional state space is displayed schematically in Fig. 7a. Here, the stable steady states A (yellow), B (blue), and C (green) are located at the minima of the landscape. Basins of attraction are displayed topographically and are the same color as their associated steady state.

In the CDI model, as demonstrated in Fig. 2 three distinct steady states can be reached from the same initial condition by administering different interventions (antibiotic administration and CD exposure). In earlier work we showed that the transition between these steady states is sudden as a function of the magnitude of the intervention (for example, an antibiotic dose of 0.70 units leads to steady state C, while a dose of 0.72 units leads to steady state E) jones2018silico . Thus, the basins of attraction of these three steady states are touching, but the structure of this high-dimensional dynamical landscape is difficult to visualize. To rectify this, attractor networks can be used to compress information about the basins of attraction of the high-dimensional phase space into a visually-digestible form. Attractor networks were originally introduced by Wang et al. to study the controllability of gene regulatory networks associated with cancer WangLai2016 .

In the schematic Fig. 7, panel (B) displays the attractor network for the dynamical landscape in panel (A). In the attractor network, nodes are steady states of the high-dimensional system, edges are convex combinations of pairs of steady states, and the colors along edges correspond to the basins of attraction along each convex combination.

Attractor networks are especially valuable for mapping the landscape near a few steady states in a high-dimensional system. For example, Fig. 8 presents an attractor network for the five reachable steady states of the CDI model (described in Fig. 3). The attractor network preserves geometric information about the basins of attraction of these five steady states, without requiring awkward visualization of an 11-dimensional state space.

The realization of antibiotic-induced CDI, as demonstrated in Fig. 2, indicated that (i) antibiotic administration shifted a microbial trajectory from the healthy steady state C towards the antibiotic-depleted steady state E, (ii) antibiotic exposure coupled with CD inoculation caused the microbial trajectory to flow towards the CD-infected steady state D, and (iii) trajectories in the basin of attraction of the CD-infected steady state D could be driven to the healthy steady state C through the introduction of a foreign population. The attractor network provides a compressed description of how an ecological state responds to external interventions, and complements and strengthens the numerical proof-of-concept analysis of antibiotic-induced CDI in Fig. 2.

The attractor network acts as a “roadmap” for a dynamical landscape, places microbial trajectories in the context of other nearby attractors, and can guide the administration of interventions in a gLV model. In discussing how to navigate gLV systems, we will assume that we are able to drive the state of the system along convex combinations of pairs of steady states.

Explicitly, represent steady states A, B, and C in Fig. 7 by the vectors $\vec{x}_{A}$, $\vec{x}_{B}$, and $\vec{x}_{C}$, respectively. Then, starting from steady state $i$, a direct intervention in the direction of steady states $j$ corresponds to a transplant $\tilde{\vec{v}}=p(\vec{x}_{j}-\vec{x}_{i})$ in the gLV equations Eq. (2), where $p$ varies between 0 and 1 and describes the magnitude of the intervention. For example, in the schematic attractor network Fig. 7b, starting from steady state A an introduced transplant $\tilde{\vec{v}}_{1}=0.3(\vec{x}_{B}-\vec{x}_{A})$ will drive the state of the system into the basin of attraction of steady state B.

Explicitly, we seek to minimize the size of the intervention required to drive a system from an initial steady state into the basin of attraction of a target steady state $\lVert\tilde{\vec{v}}\rVert$. Using the same schematic attractor network from Fig. 7b, driving an initial condition at steady state B into the basin of attraction of steady state C will require a transplant $\tilde{\vec{v}}_{2}=0.3(\vec{x}_{C}-\vec{x}_{B})$. In this schematic the distances between each pair of steady states is defined to be one unit, so applying these two interventions in a sequential manner will drive an initial condition at steady state A into the basin of attraction of steady state C with a total intervention size of 0.6. Alternatively, to drive an initial condition at steady state A into the basin of attraction of steady state C directly requires a transplant $\tilde{\vec{v}}_{3}=0.7(\vec{x}_{C}-\vec{x}_{A})$, with a total intervention size of 0.7. Thus, constructing an attractor network demonstrates that inherent ecological dynamics can be leveraged to efficiently control the state of an ecological system.

In the CDI model, the attractor network in Fig. 8 reveals that scenarios exist in which the most efficient control method follows indirect paths. In Fig. 9 the microbial time courses associated with two control protocols— one that follows a direct path, and another that follows a circuitous path— are compared.

Explicitly, let steady state $i$ of the CDI model correspond to the vector $\vec{y}_{i}$. To transition from the antibiotic-depleted steady state E directly to the healthy steady state C requires a transplant $\vec{v}_{1}=0.162(\vec{y}_{C}-\vec{y}_{E})$, with total transplant size $\lVert v_{3}\rVert=2.576$, where $\lVert\cdot\rVert$ is the 2-norm and each transplant is in units of $10^{11}$ microbes. A similar control protocol is plotted in Fig. 9a, which successfully drives the state of the system into the basin of attraction of steady state C (the plotted transplant $\vec{w}_{1}$ is chosen to be slightly larger than the required transplant $\vec{v}_{1}$ for demonstrative purposes).

However, in this CDI model it is most efficient to apply two sequential interventions $\vec{v}_{2}$ and $\vec{v}_{3}$: the first transplant drives the system towards the CD-infected steady state D, $\vec{v}_{2}=0.001(\vec{y}_{D}-\vec{y}_{E})$; and the second transplant drives the system towards the healthy steady state C, $\vec{v}_{3}=0.022(\vec{y}_{C}-\vec{y}_{D})$. Here, since steady state E is unstable in the direction of steady state D, an infinitesimal transplant is all that is needed (though for practical purposes we use a value of 0.001). The sizes of these two sequential transplants are $\lVert v_{2}\rVert=0.0001$ and $\lVert v_{3}\rVert=0.360$. A similar sequential control protocol is demonstrated in Fig. 9b; once again, for demonstrative purposes the transplants $\vec{w}_{2}$ and $\vec{w}_{3}$ are chosen to be slightly larger than the required transplants $\vec{v}_{2}$ and $\vec{v}_{3}$.

Taken together, the circuitous control protocol E $\to$ D $\to$ C (total transplant size $0.36$) requires a smaller total intervention than the direct path E $\to$ C (total transplant size $2.57$). In general, taking advantage of this “ecological inertia” may reduce the magnitude of the intervention required to drive the system towards a target state.

Finally, attractor networks can be approximated in constant time (in algorithmic complexity terms) with the dimensionality-reduction technique SSR. The attractor network for the CDI model in Fig. 8 was generated numerically by simulating a set of microbial trajectories, each originating at initial conditions along convex combinations of steady-state pairs, and then tracking the steady-state outcome of each simulation. This procedure is computationally expensive even when bisection-type algorithms are implemented to identify the location of the separatrix.

As derived in previous work, SSR allows for the separatrix of the approximate SSR-reduced 2D gLV system to be computed analytically in a Taylor series about a semi-stable fixed point $(x_{a}^{*},x_{b}^{*})$ of the reduced system JonesCarlson2019 . For example, the separatrix of the 2D gLV system in Fig. 5 was created in this fashion. Though the current attractor network is relatively small, consisting of five nodes and ${5\choose 2}$ edges, scaling up the size of a numerically-calculated attractor network quickly becomes computationally infeasible. On the other and, larger attractor networks can be efficiently approximated with SSR: as a demonstration of the accuracy with which SSR captures the location of the separatrices, in Fig. 10 the squares on each edge represent the separatrix predicted by applying SSR to the corresponding steady-state pair.

The SSR-generated attractor network incorrectly predicts that convex combinations of steady states B and D will flow towards steady state D, when in fact they will flow towards steady state A. This occurs because the pairwise connections in the graph representation of the attractor network may oversimplify the complex topography of the original landscape. Still, the strong agreement between the exact and SSR-generated separatrices indicates the potential of SSR to efficiently generate large attractor networks to first-order that intuitively describe relevant ecological dynamics.

Microbiome research has been the recent beneficiary of advances in experimental 16S pyrosequencing techniques that have revealed similarities in the microbiome compositions of people suffering various diseases, which differ from the microbiome compositions of healthy individuals shreiner2015gut . These disease-associated microbiome compositions are called dysbiotic, and are observed in people that suffer cardiovascular disease, ulcerative colitis, and irritable bowel disease lee2010has ; shreiner2015gut .

In this chapter we studied a mathematical model of C. difficile infection (CDI). CDI occurs when the bacteria C. difficile (CD) colonizes the gut and becomes sufficiently abundant to induce colon inflammation and diarrhea via the production of toxins TcdA and TcdB voth2005clostridium . This disease is an explicit example of how impaired microbiome compositions cause disease.

Traditional treatments for CDI seek to eliminate the presence of CD by administering antibiotics like vancomycin or metronidazole mcfarland2002breaking . However, CDI recurrence rates that range from 30-65% and concerns about antibiotic resistance have heralded fecal microbiota transplantation (FMT) as an alternative treatment for CDI brandt2012long ; spigaglia2011multidrug . FMT is a bacteriotherapy that attempts to alter the composition of a dysbiotic microbiome (e.g., the microbiome of a person suffering CDI) by engrafting a foreign microbial population provided by a healthy donor into it. When successful, FMT causes the restoration of a healthy microbiome, which suppresses the growth of CD as well as the production of CDI-associated toxins khoruts2016understanding . The success rate of FMT for treating CDI approaches 90% brandt2012long .

FMT has been proposed as a treatment for other gastointestinal diseases including irritable bowel syndrome (IBS), inflammatory bowel disease, ulcerative colitis, and Crohn’s disease. Recent clinical trials have returned promising but inconclusive evidence regarding the success of FMT in treating these conditions coryell2018gut ; franzosa2018gut ; de2013transient ; anderson2012systematic . For other conditions with distinct microbial signatures— for example, autism spectrum disorder (ASD) or metabolic syndrome vuong2017emerging ; perry2016acetate — a more intrinsic question regarding the efficacy of bacteriotherapies remains: can symptoms of these conditions be ameliorated by altering the gut microbiome composition? For ASD, preliminary evidence supports this hypothesis. Patients with ASD were given Microbiota Transfer Therapy, a treatment that consists of an initial course of antibiotics followed by regular FMT treatments for ten weeks, and this treatmenet caused lasting shifts to their microbiome compositions and reductions in their ASD symptoms kang2019long .

At this time, the ability of FMT to treat microbiome-associated diseases is variable. FMT appears to be an effective treatment for CDI, but recent clinical studies have been unable to conclusively show that FMT is an effective treatment for IBS xu2019efficacy . Part of this discrepancy might be due to antibiotics, which are typically administered to CDI patients before FMT is attempted. More generally, the factors that contribute to the success or failure of FMT are not yet fully known. To shed light on these factors, this chapter used mathematical models to examine the ecological dynamics that underlie FMT.

Models of the microbiome can span levels of organization, stochasticity, neutrality, and complexity. Recent advances in multi-omics sequencing have buoyed microbiome modeling by measuring the abundances of the microbial genomes, mRNA transcripts, proteins, and metabolites contained in the microbiome biggs2015metabolic .

For microbial ecosystems with well-characterized metabolic pathways, flux balance analysis (FBA) predicts the production and consumption of metabolites by microbes at steady state using inferred stochiometric matrices biggs2015metabolic . In dynamic flux balance analysis (dFBA), FBA is generalized to track metabolite abundances over time, at the expense of requiring additional kinetic parameters in the model biggs2015metabolic . In both of these metabolite-based analyses, fine-grained details about metabolite kinetics need to be either measured or fit. This required level of detail makes FBA and dFBA ill-suited for analyzing complex microbial ecosystems like the human microbiome which, for example, consists of $\sim$1000 microbial species and $\sim$100000 metabolites lloyd2016healthy ; wishart2018drugbank .

To account for the complexity present in real-world microbial systems, coarse-graining approaches have been employed to model aggregate microbial dynamcs and processes. For example, in wastewater treatment activated sludge models (ASMs) track the abundance of only a few nutrients (e.g. organic carbon, nitrogen, and phosphorus) and the population dynamics of only a few aggregate microbial populations (e.g. microbes that consume carbon, nitrogen, and phosphorus) bucci2014towards . Though less detailed than FBA and dFBA models, ASMs require significantly fewer parameters and are able to provide relevant information regarding the design and function of wastewater treatment processes.

Generalized Lotka-Volterra (gLV) models, which are studied in this chapter, assume that the production and consumption of nutrients can be implicitly captured by pairwise interactions between microbial populations. Thus, gLV models only track microbial population abundances. Some gLV models (including the CDI model) assume that microbial species within a genus are indistinguishable, and study microbial dynamics coarse-grained at the genus level. Part of the convenience of gLV models is that they may be parameterized with time-series microbial abundance data, which is readily available due to advances in genomic sequencing. Several algorithms have been developed to infer growth rate and interaction parameters of gLV models from time-series abundance data (e.g. LIMITS fisher2014identifying ; bucci2016mdsine ). gLV models have been used to investigate microbial interactions in cheese mounier2008microbial , to probe how antibiotic perturbations alter the gut microbiome in the context of CDI stein2013ecological , and to inform treatment of polymicrobial urinary tract infections de2017interaction . When gLV models accurately approximate the underlying microbial systems, simulated interventions (e.g. antibiotic administration or fecal microbiota transplantation) can guide clinical efforts to alter the composition of diseased microbiomes stein2013ecological ; buffie2015precision ; jones2018silico .

Each of these previously discussed models make the simplifying assumption that microbial processes and the resulting microbial dynamics are deterministic. In ecology these deterministic approaches are typically used to model niche processes. While microbial dynamics can be relatively consistent in some systems (e.g. in soil microbial communities friedman2017community ), variability— reflecting neutral processes— is generally evident in microbial communities (e.g. in human li2016testing or fly obadia2017probabilistic microbiomes). Mathematical models have been developed to account for this variability faust2018signatures . For example, Sloan et al. proposed a model to predict microbial abundances based solely on birth-death processes and immigration following Hubbell’s neutral theory sloan2006quantifying , while the Ricker model has been used to introduce stochasticity into gLV models fisher2014identifying .

Existing theory describes the controllability of polynomial differential equations using differential geometry Kunita1979 , and these techniques should be applicable to gLV systems. For example, the Pontryagin maximum principle provides a necessary condition for whether a proposed control protocol is optimal BoltyanskiyPontryagin1962 . This approach has been previously applied to a optimal sampled-data control problem in the context of a biomechanical model that predicts how muscles respond to electrical stimulations BourdinTrelat2016 . In this biomechanical study, a version of the Pontryagin maximum principle specific to optimal sampled-data control problems was used to identify the optimal timing of electrical stimulations that maximize the muscular force response BakirRouot2020 . Since our microbiome control problem also treats direct interventions as Dirac delta inputs, this modified version of the Pontryagin maximum principle could be applied to gLV systems.