# A mathematical model of B. thetaiotaomicron, M. smithii, and E. rectale
interactions in the human gut

Melissa A. Adrian  Bruce P. Ayati  University of Iowa, Department of
Mathematics, Iowa City, IA, 52242, USA  Ashutosh K. Mangalam  University of
Iowa, Department of Pathology, Iowa City, IA, 52242, USA

###### Abstract

The human gut microbiota is a complex ecosystem that affects a range of human
physiology. In order to explore the dynamics of the human gut microbiota, we
used a system of ordinary differential equations to mathematically model the
biomass of three microorganism populations: Bacteroides thetaiotaomicron,
Eubacterium rectale, and Methanobrevibacter smithii. Additionally, we modeled
the concentrations of relevant nutrients necessary to sustain these
populations over time. Our model highlights the interactions and the
competition among these three species. These three microorganisms were
specifically chosen due to the system’s end product, butyrate, which is a
short chain fatty acid that aids in developing and maintaining the intestinal
barrier in the human gut. The basis of our mathematical model assumes the gut
is structured somewhat similar to a chemostat, with bacteria and nutrients
exiting the gut at a rate proportional to the volume of the chemostat, the
rate of volumetric flow, and the biomass or concentration of the particular
population or nutrient. We performed global sensitivity analyses using Sobol’
sensitivities to estimate the relative importance of model parameters on
simulation results.

## Introduction

The human gut microbiota is the collection of microorganisms located in the
stomach, large intestines, and small intestines, and this system plays an
important role in sustaining overall human health. Each individual’s gut
composition is unique, and, besides the host’s genetic makeup, their long-term
dietary patterns, specifically concerning the types and amounts of
carbohydrates, proteins, and fats consumed, affect its composition [1, 2].
Specifically, gut microbiota have been shown to be in involved in numerous
physiological processes, such as digestion of undigested food (complex
starch), development and regulation of immune system, blocking growth of
pathogens, and generation of neurotransmitter and vitamins. Thus, any changes
in environmental conditions in this gut ecosystem (microbiota) can result in a
shift in its composition, which can predispose and/or worsen chronic
inflammatory diseases such as inflammatory bowel diseases, multiple sclerosis,
rheumatoid arthritis and neurological diseases like Alzheimer’s disease,
Parkinson disease and autism [3, 4, 5]. Though the changes in the bacterial
composition of the human gut microbiota have been shown to be associated with
rapid changes in metabolism and overall health, the underlying interactions
among species within this microbiota are not yet well understood [2]. Among
all the factors linked with regulation of gut microbiota, diet has emerged as
the strongest as it can override genetic influences on microbiota [6, 7]. A
better understanding of the microbial dynamics, particularly mathematical
approaches, may elucidate the role that the gut microbiota plays in human
diseases [1].

Common approaches to analyzing the human gut microbiota for its diversity and
composition include, but are not limited to, next-generation sequencing,
metatranscriptomics, culturomics, and mass spectrometry analyses [8, 2]. These
approaches give little insight as to how species interact amongst themselves
and with the human host [2]. Mathematical modeling, however, attempts to
answer such questions and can supplement the knowledge gained from these
common approaches. Mathematical models can provide evidence to strengthen a
hypothesis about these interactions and offer guiding principles for further
study of a phenomenon [9].

We use ordinary differential equation (ODE)-based modeling as the main tool
for our analysis, which tracks information about biomass and concentration
levels over time rather than genomic information as in the case of GEMs. In
this so-called dynamical system approach, we can identify the system’s driving
parameters and analyze its stability. Despite the advantage of having easily
interpretable terms in the model’s equations, this approach still has some key
drawbacks, specifically the determination of unknown parameters.

The parameters used in biological models like ours include substrate
conversion rates, utilization rates, and cellular death rates, which typically
can be determined and validated by laboratory experiments. However, some gut
microorganisms are unable to be cultivated in a laboratory setting, leaving
specific information about these species’ metabolic activity unknown.
Additionally, ODE models are often used to track only a small subset of
species within a community. While this is often due to the difficulty of
determining a larger number of parameters when the system is expanded [8],
focusing on a few subspecies at a time may aid in identifying the dominant
relationships in that subsystem.

Despite the impracticality of ODE-based modeling for some very large systems,
ODE models paired with other types of modeling, such as agent-based modeling,
can provide an insightful understanding of the gut microbiota’s dynamics and
interactions [8]. As a first step in our modeling efforts with its limitations
in mind, we consider an ODE-based modeling approach for a small-scale system
of three abundant microorganisms in the human gut microbiota. Our ODE model is
structured similar to models of chemostats, which are idealized laboratory
systems for microbial ecology that have a long history of mathematical
modeling behind it [10]. One contribution of this paper is the use of a
similar approach, rather than the “batch” models in other gut microbiome
models [8], so that the vital role of the inflows and outflows in the human
gut ecosystem has a more faithful representation.

We focus our representation on a subsystem of the human gut microbiota. The
three microorganisms Bacteroides thetaiotaomicron, Eubacterium rectale, and
Methanobrevibacter smithii were chosen to be the focus of our model. These
species represent the three main phyla in the human gut: Bacteroidetes,
accounting for 17-60% of the total biomass; Firmicutes, 35-80% of the biomass;
and Euryarchaeota the bulk of the remainder[11]. The system’s main product,
butyrate, is of specific interest due to its role in sustaining human health.
Butyrate provides energy to colonocytes, affects overall energy homeostasis,
and inhibits histone deacetylase, which is an enzyme that directly affects
colorectal cancer[11]. Along with butyrate, other notable intermediates and
products in this system include acetate, propionate, glutamine, carbon dioxide
(CO2), hydrogen (H2), and methane (CH4). Acetate, propionate, and butyrate,
which are short chain fatty acids, are absorbed in the gut’s epithelial cells
and regulate an individual’s immune system and metabolism[11]. Individually,
acetate acts as a substrate for cholesterol synthesis and lipogenesis [11];
propionate regulates gluconeogenesis and cholesterol synthesis[11]; glutamine
fuels the metabolism and maintains the intestinal barrier[12]; and the gases
CO2, H2, and CH4 are products of bacterial fermentation that can cause
intestinal discomfort in excess[13]. As for the microorganisms themselves, M.
smithii removes hydrogen gas, which affects bacterial fermentation and energy
gathering, and produces methane gas[11]; E. rectale produces butyrate, which
is beneficial to the gut’s epithelial cells[11]; and B. thetaiotaomicron
utilizes dietary polysaccharides and indirectly facilitates butyrate
production with its outputs [2].

In our model, we represent the interactions of three short chain fatty acid
producing/utilizing gut microorganisms (see Figure 1) and their metabolites.
B. thetaiotaomicron is an abundant bacterial species in the human gut
microbiome whose main function is the utilization of polysaccharides [14, 2].
Through polysaccharide degradation, B. thetaiotaomicron contributes to the
overall ecosystem diversity in the colon, which is its regular environment
[15]. B. thetaiotaomicron can survive solely on the uptake of carbon-rich
polysaccharides[16]; however, its growth is enhanced in the presence of
inorganic ammonia due to inorganic ammonia’s contribution of nitrogen to B.
thetaiotaomicron’s metabolism [17]. Through the utilization of inorganic
ammonia, B. thetaiotaomicron can synthesize all amino acids that are essential
to human health [2], which makes this bacteria a focus of study.

When both E. rectale and B. thetaiotaomicron are present in an environment, B.
thetaiotaomicron up-regulates gene expression for starch utilization and the
degradation of specific glycans that E. rectale is unable to utilize.
Simultaneously, E. rectale down-regulates the genes associated with glycan
degradation even though it cannot grow efficiently without a carbohydrate
source. Previous research on the interactions of these two species suggests
that the presence of B. thetaiotaomicron enhances E. rectale’s ability to
uptake nutrients[18]. E. rectale shifts from uptaking polysaccharides to
utilizing amino acids, such as glutamine, when B. thetaiotaomicron is
present[2].

M. smithii, which is one of the main methanogenic archaeon in the human gut,
improves the productivity of carbohydrate metabolism by utilizing H2 from E.
rectale and formate from B. thetaiotaomicron to produce methane gas. This
process prevents the environment from becoming too saturated with B.
thetaiotaomicon and E. rectale’s by-products, which consequently improves
carbohydrate metabolism. Additionally, M. smithii removing H2 in this
environment allows for B. thetaiotaomicron to generate NAD+, which is used for
glycolysis, a fundamental process in producing cellular energy [18].

Figure 1: Graphical representation of the interactions of the three species,
B. thetaiotaomicron, E. rectale, and M. smithii, and their metabolites. Each
microorganism or substance given a letter is included in the mathematical
model equations. This figure was developed from information presented in Ji &
Nielsen 2015 [2], Shoaie et al. 2013[11], and Adamberg et al. 2014 [15].

The information known about these three species’ interactions is translated
into a graphical representation in Figure 1. This schematic highlights the
competition among these species, specifically between E. rectale and M.
smithii and between B. thetaiotaomicron and E. rectale.

The results of our modeling effort suggest the efficacy of a relatively simple
ODE approach, similar to those used for chemostats, in understanding and
potentially predicting the dynamics of important subsystems in the larger
human gut ecosystem.

## Results

Using our mathematical model of Bacteroides thetaiotaomicron, Eubacterium
rectale, and Methanobrevibacter smithii and their metabolites in equations 1
and 2, we present numerical results, estimation of key parameters, and
sensitivity analysis.

### Numerical Results

Using the set of parameter estimates in Table 5, the solutions to our full
model are given for the dynamics of B. thetaiotaomicron, M. smithii, and E.
rectale in Figure 2, for acetate in Figure 3, for CO2 and H2 in Figure 4, and
for polysaccharides in Figure 5. These plots show the solutions to our full
system of ODEs in equations (1) and (2) from 1,000 to 1,100 hours. This time
range was selected in order to allow for a sufficient amount of time to pass
in order for the system to converge to a regular oscillation. The center
values of the oscillations are represented by dotted lines in each plot. In
Figures 3, 4, and 5, the data values given in Table 3 are superimposed on
their respective plots as dashed lines in order to provide a visual
representation of the error between the mean of the regular oscillation and
the observed data. Despite some error in the center value of the ODE solutions
compared to the data, the observed data values for acetate, polysaccharides,
and the gases CO2 and H2 are contained in the oscillation range of the ODEs’
numerical solutions, so we conclude that our parameter estimates sufficiently
fit the data. Thus, the relatively simple model we have presented here is
consistent with the data in our literature review.

Figure 2: Plot of the solutions to equations (1a), (1b), and (1c) using the
observed data in Table 3 over the time interval 1,000 to 1,100 hours. Based on
Table 3, the sum of the three species’ biomass should be 0.001412 gDW. The
middle of the solutions to the ODEs using the parameter estimates in Table 5
for all three species sums to 0.000925 gDW, resulting in an error of
$-$0.000487 gDW. Figure 3: Plot of the solutions to equation (2a) using the
observed data in Table 3 over the time interval 1,000 to 1,100 hours. Based on
Table 3, the center of acetate’s oscillations should be 9.71 $\mu$M. The
middle concentration of the ODE solutions for acetate is 10.47 $\mu$M,
resulting in a 0.76 $\mu$M error. Figure 4: Plot of the solutions to equation
(2b) using the observed data in Table 3 over the time interval 1,000 to 1,100
hours. Based on Table 3, the center of CO2 and H2’s oscillations should be
roughly 7.96 $\mu$M. The middle concentration of the ODE solutions for CO2 and
H2 is 7.80 $\mu$M, resulting in a $-$0.16 $\mu$M error. Figure 5: Plot of the
solutions to equation (2c) using the observed data in Table 3 over the time
interval 1,000 to 1,100 hours. Based on Table 3, the center of
polysaccharides’ oscillations should be roughly 32.06 $\mu$M. The middle
concentration of the ODE solutions for polysaccharides is 29.50 $\mu$M,
resulting in a $-$2.56 $\mu$M error.

### Sensitivity Analysis

Based on the fitted parameter values and the available data from the
literature, we analyzed our model’s sensitivity to its parameters. In order to
identify the parameters in the model that have the greatest effect on the
model output, we conducted sensitivity analysis on our model by computing the
first- and total-order effects based on Sobol’ indices.

Sensitivity analysis can be defined as the study of how uncertainty in the
model input propagates into uncertainty in the model output[19]. Once
parameters are estimated, sensitivity analysis can be used to identify the
driving parameters of the system that contribute to the most variability in
the model’s output. Sensitivity analysis tests the robustness of the model,
identifies if the model relies on weak assumptions, and allows for model
simplification[19].

The results of our sensitivity analysis using the methods described in the
section titled "Sensitivity Analysis Methods" are given in Tables 7 and 7.

| $k$ | $\theta$ | Mean | Variance
---|---|---|---|---
$\beta_{a}$ |  $5.00$$\times 10^{6}$ | $2.00$ |  $1.00$$\times 10^{7}$ |  $2.00$$\times 10^{7}$
$\beta_{B}$ | $0.60$ | $2.00$ | $1.20$ | $2.40$
$\beta_{E_{1}}$ | $0.40$ | $2.00$ | $0.80$ | $1.60$
$\beta_{E_{2}}$ | $0.30$ | $2.00$ | $0.60$ | $1.20$
$\beta_{h_{1}}$ | $75.00$ | $2.00$ | $150.00$ | $300.00$
$\beta_{h_{2}}$ |  $1.67$$\times 10^{4}$ | $2.00$ |  $3.33$$\times 10^{4}$ |  $6.67$$\times 10^{4}$
$\beta_{M_{1}}$ | $0.40$ | $2.00$ | $0.80$ | $1.60$
$\beta_{M_{2}}$ | $0.25$ | $2.00$ | $0.50$ | $1.00$
$\beta_{p}$ |  $5.00$$\times 10^{3}$ | $2.00$ |  $1.00$$\times 10^{4}$ |  $2.00$$\times 10^{4}$
$\gamma_{a}$ | $100.00$ | $2.00$ | $200.00$ | $400.00$
$\gamma_{B}$ | $5.00$ | $2.00$ | $10.00$ | $20.00$
$\gamma_{h}$ | $75.00$ | $2.00$ | $150.00$ | $300.00$
$\gamma_{p}$ | $200.00$ | $2.00$ | $400.00$ | $800.00$
$\mu_{aE}$ |  $1.25$$\times 10^{3}$ | $2.00$ |  $2.50$$\times 10^{4}$ |  $5.00$$\times 10^{4}$
$\mu_{aM}$ |  $2.50$$\times 10^{4}$ | $2.00$ |  $5.00$$\times 10^{4}$ |  $1.00$$\times 10^{5}$
$\mu_{hM}$ | $20.00$ | $2.00$ | $40.00$ | $80.00$
$\mu_{pB}$ |  $1.00$$\times 10^{5}$ | $2.00$ |  $2.00$$\times 10^{5}$ |  $4.00$$\times 10^{5}$
$\mu_{pE}$ |  $2.50$$\times 10^{3}$ | $2.00$ |  $5.00$$\times 10^{3}$ |  $1.00$$\times 10^{4}$
| $\alpha$ | $\beta$ | Mean | Variance
$q$ | $5.00$ | $95.00$ | $0.05$ |  $4.70$$\times 10^{-2}$
Table 1: Table of specified distributions for the parameters in our full
model, including the mean and variance of each distribution. The parameters
whose distribution is specified with $k$ and $\theta$ follow a
Gamma($k,\theta$) distribution. The parameter $q$, however, follows a
Beta($\alpha,\beta$) distribution.

Using our code in Appendix A to produce the first- and total-order Sobol’
indices for our full ODE model, we obtained the results seen in Tables 7 and
7. Based on Table 7 of the first-order Sobol’ indices, $\beta_{B}$ is an
influential parameter for B. thetaiotaomicron’s output variance;
$\beta_{E_{2}}$ for E. rectale’s variance; $\beta_{M_{2}}$ for M. smithii’s
variance; no first-order parameter for acetate’s variance; $\beta_{B}$ and
$\beta_{M_{2}}$ for CO2 and H2’s variance; and $\beta_{B}$ for
polysaccharides’ variance. Based on Table 7 of the total-order Sobol’ indices,
$\beta_{B}$ is an influential parameter for B. thetaiotaomicron’s output
variance; $\beta_{E_{2}}$ and $\beta_{B}$ for E. rectale’s variance;
$\beta_{M_{2}}$, $q$, and $\beta_{B}$ for M. smithii’s variance; $\beta_{B}$,
$\beta_{E_{1}}$, $\beta_{E_{2}}$, $\beta_{M_{1}}$, $\beta_{M_{2}}$, and $q$
for acetate’s variance; $\beta_{B}$, $\beta_{E_{2}}$, $\beta_{M_{1}}$,
$\beta_{M_{2}}$, and $q$ for CO2 and H2’s variance; and $\beta_{B}$,
$\beta_{E_{2}}$, and $q$ for polysaccharides’ variance. Based on the changes
in the estimated indices from first-order to total-order indices, there
appears to be a higher-order interaction among model parameters for E.
rectale’s, M. smithii’s, acetate’s, CO2 and H2’s, and polysaccharides’ output
variance. We note that observing Sobol’ indices whose sum is greater than 1,
as seen in Table 7, is evidence of potentially correlated inputs in the model
[20]. Overall, the two main driving parameters of the ODE model appear to be
the reproduction coefficient for B. thetaiotaomicron ($\beta_{B}$) and the
ratio of the volume of the chemostat over the volumetric flow ($q$). Because
these parameters are identified to be sensitive to perturbations, careful
attention needs to be given to estimating these parameters and reducing the
error of their estimates in future efforts.

## Discussion

We suggest and illustrate that a mathematical representation similar to those
for chemostats is a natural way to capture the effects of the inflow and
outflow in the gut on its microbial dynamics.

The three species chosen to be the focus of our model, B. thetaiotaomicron, E.
rectale, and M. smithii, play an important role in polysaccharide degradation
and the production of buyrate, which both aid in the human gut’s ability to
absorb nutrients through the epithelial cells[11]. The system of the three
microbial species has been considered in previous works[11, 2], which largely
informed our knowledge of this system. By creating a mathematical model based
on the interactions of these species, we have analyzed their interactions and
identified aspects of this system that should be further explored through
empirical investigation.

Due to the limited availability of data for a more rigorous parameter
estimation, our sensitivity analysis took on additional importance. Through
the results of the sensitivity analysis using first- and total-order Sobol’
indices, we more narrowly identified specific links in the microbial food web
that would be fruitful targets for additional empirical work. Specifically, we
identified the growth rate coefficients $\beta_{B}$, $\beta_{E_{1}}$,
$\beta_{E_{2}}$, $\beta_{M_{1}}$, $\beta_{M_{2}}$, and the scaled volumetric
flow ($q$) as being largely significant in contributing to the variance of the
model output, including higher-order interactions among these parameters. With
these results, we suggest that estimates of these significant parameters be
obtained through laboratory experimentation in order to capture these values
to a higher degree of precision and accuracy.

The significant growth rate coefficients ($\beta$ parameters) correspond to
the growth rates of the three microorganisms when supported on a medium
containing specific nutrients. Experiments should focus on cultivating these
microorganisms in isolation in germ-free mice with only one nutrient [11].
Despite the fact that these microorganisms are able to be cultured in
isolation of other microorganisms, they may not be able to be sustained on one
single nutrient, but rather require the presence of additional substrates. In
this case, nonlinear effects from these secondary nutrients would factor in to
the resulting estimates of the growth parameters. Table 2 provides a general
outline of the the specific nutrient and microorganism necessary to estimate
each parameter. For example, the growth rate coefficient for B.
thetaiotaomicron ($\beta_{B}$) can be experimentally estimated by cultivating
B. thetaiotaomicron in a medium containing only polysaccharides. Additionally,
experiments approximating the rate of flow of the digestive system’s fluids
would largely inform the estimate for the true value of the volumetric flow
rate parameter $q$.

Parameter | Microorganism | Nutrient
---|---|---
$\beta_{B}$ | B. thetaiotaomicron | polysaccharides
$\beta_{E_{1}}$ | E. rectale | acetate
$\beta_{E_{2}}$ | E. rectale | polysaccharides
$\beta_{M_{1}}$ | M. smithii | acetate
$\beta_{M_{2}}$ | M. smithii | CO2 and H2
Table 2: Table of microorganism growth rate parameters suggested to be
estimated experimentally.

Due to the possibility of correlation among model parameters, variance-based
sensitivity analyses specifically for correlated parameters should be explored
and applied to this system, such as the methods discussed in Iooss and Prieur
2019[20] and Rabitz 2010[21]. In regards to specifying the prior distributions
on parameters, further analyses should include testing differing prior
parameter distributions in order to compute the Sobol’ indices, especially to
explore the amount of variability in results based on the chosen parameter
distributions[19].

Our main suggestion for future work is to collect more complete longitudinal
data on this biological system, including for all three species B.
thetaiotaomicron, E. rectale, and M. smithii and its relevant substrates. With
this experimental data, more precise estimates of the model parameters can be
achieved. With this approach, particular attention can be paid to estimating
the parameters that were identified as sensitive by the Sobol’ indices.

Additionally, further extensions of our model may include relaxing some of the
assumptions of a general chemostat model. Specifically, the rate of volumetric
flow through the chemostat can be generalized to reflect aspects of the
natural flow of the gut, such as periodically restricted flow, and the
framework of a single vessel representation can be extended to consider
additional compartments to mimic the system of digestive organs in the human
body. Aside from the assumptions of a simple chemostat model, future
implementations of a similar model can account for absorption rate of
substrates within the gut, which is especially important to consider for
substances like amino acids. With these potential improvements of our baseline
model, additional aspects about the dynamics of this biological system can be
uncovered, and these improvements to our model could fuel further research
directions related to this system.

Our model has relevance to human microbiome studies as it can be used to test
the significance of clinical findings. For example, microbiome studies for
human autoimmune disease, such as multiple sclerosis, have shown modulation of
butyrate and methane-producing gut bacteria [22]. This can be incorporated
into our modeling framework in combination with experimental studies to
determine the mechanism through which gut bacterial modulate disease.

## Methods

The key assumption underlying our mathematical model is that the human gut can
be represented as a chemostat. A chemostat is a laboratory device used in the
simulation and ecological study of populations, which provides an idealized
representation of naturally occurring phenomena and has a rich history of
mathematical representation. Though the conditions of a chemostat are
simplified and controlled in a laboratory setting, a chemostat can be useful
in the study of population dynamics and the underlying mechanisms of
interactions among populations. The simple chemostat model is a first step in
developing an initial theoretical framework, which can then be refined and
extended[10].

In developing an ODE model to represent our chosen microbiota subsystem, we
assume that the human gut acts in a manner similar to that of a chemostat,
which is an approach not used, to our knowledge, in other work. We utilize
ODE-based dynamical systems modeling to track the changes in biomass of the
three prevalent microorganisms, B. thetaiotaomicron, E. rectale, and M.
smithii, as well as their chemical inputs, intermediates, and byproducts, with
the goal of providing a better understanding of their interactions within this
subsystem.

In order to supplement the knowledge gained from transcriptomic analysis and
GEMs, we incorporate this information into our own mathematical interpretation
of this small-scale system using ODEs. Throughout our model explanation, we
present information learned through previous investigations of B.
thetaiotaomicron, E. rectale, and M. smithii that we utilized in creating our
mathematical model.

In constructing a first chemostat model, we assume that the contents of the
vessel are well-mixed, the rate at which liquid enters the system equals the
rate at which the well-stirred contents leave the compartment, and that some
significant factors potentially affecting growth, such as temperature, are
held constant[10].

For specificity we assume that microorganisms grow at a rate following the
Monod form

$\beta_{X}\left(\frac{u}{u+\gamma}\right)X,$

where $\beta_{X}$ is the maximum birth rate of population $X$, $u$ is the
concentration of the nutrient population on which $X$’s growth depends,
$\gamma$ is the Michaelis-Menten constant, and $X$ is the concentration of the
microorganism[23]. This form is commonly used, especially for a first effort.
Other functional forms that could be used include Hill functions[24], which
are generalizations of the Monod form, and S-forms[25]. The accompanying large
increase in the number of parameters, and lack of data to capture the detail
these parameters provide, means that they are less appropriate for our effort
here than the Monod form.

### Model Variables

The variables in our model for the microbial and chemical species depend only
on one independent variable, time, denoted by $t$. These dependent variables
are

  * •

$B$ for B. thetaiotaomicron,

  * •

$E$ for E. rectale,

  * •

$M$ for M. smithii,

  * •

$a$ for acetate,

  * •

$h$ for CO2 and H2,

  * •

$p$ for polysaccharides.

### Model Equations

The equations for the three microbial species are

$\displaystyle\frac{dB}{dt}$
$\displaystyle=\beta_{B}\Psi_{\gamma_{p}}(p)B-qB,$ (1a)
$\displaystyle\frac{dE}{dt}$
$\displaystyle=\left[\beta_{E_{1}}\Psi_{\gamma_{a}}(a)+\beta_{E_{2}}\left(1-\Psi_{\gamma_{B}}(B)\right)\Psi_{\gamma_{p}}(p)\right]E-qE,$
(1b) $\displaystyle\frac{dM}{dt}$
$\displaystyle=\left[\beta_{M_{1}}\Psi_{\gamma_{a}}(a)+\beta_{M_{2}}\Psi_{\gamma_{h}}(h)\right]M-qM.$
(1c)

The equations for the three chemical species are

$\displaystyle\frac{da}{dt}$
$\displaystyle=\beta_{a}\Psi_{\gamma_{p}}(p)B-qa-[\mu_{a,E}E+\mu_{a,M}M]a,$
(2a) $\displaystyle\frac{dh}{dt}$
$\displaystyle=\beta_{h_{1}}\Psi_{\gamma_{a}}(a)E+\beta_{h_{2}}\Psi_{\gamma_{p}}(p)B-qh-\mu_{h,M}hM,$
(2b) $\displaystyle\frac{dp}{dt}$
$\displaystyle=\beta_{p}q(\textrm{cos}(t)+1)^{3}-qp-[\mu_{p,B}B+\mu_{p,E}E]p.$
(2c)

The $\Psi$ and $q$ terms in these equations are defined by

$\displaystyle\Psi_{\gamma}(u)$ $\displaystyle=\frac{u}{u+\gamma},$ (3)
$\displaystyle q$ $\displaystyle=\frac{V}{Q}.$ (4)

In our model we assume that there is a high enough rate of turnover of fluids
in the human gut such that these microorganisms are almost always flushed out
of the system before their life expectancy, so death terms are neglected in
these three equations. All three microorganism equations thus have the general
format:

$\Delta_{i}=P_{ij}-F_{i},$

where $\Delta_{i}$ is the rate of change of biomass for microorganism $i$,
$P_{ij}$ is the rate at which microorganism $i$ proliferates based on the
availability of substance $j$, and $F_{i}$ is the rate at which microorganism
$i$ is flushed out of the gut. In all three microorganism equations, the term
$F_{i}$ is the biomass of the given population multiplied by the constant $q$.
The fixed quantity $q$ is interpreted as the rate at which the contents of the
gut leave the system as expressed in equation (4), where $V$ is the volume of
the chemostat and $Q$ is the rate of volumetric flow within the chemostat[10].

In order for E. rectale to grow in biomass, acetate or polysaccharides need to
be present in the ecosystem[11]. We assume that E. rectale has different
maximum growth rates depending on each nutrient, leading us to split
$\beta_{E}$ into three different, related constants. Because E. rectale shifts
to uptaking inorganic ammonia when B. thetaiotaomicron is present, we included
the term $\big{(}1-\Psi_{\gamma_{B}}(B)\big{)}$ to reflect this shift. M.
smithii depends on the presence of acetate and the gases CO2 and H2, which are
both incorporated in the standard Monod form. Again, we assume that M. smithii
grows at different maximum growth rates in the presence of only one of these
metabolites, which lead to the separation of $\beta_{M}$ into the constants
$\beta_{M_{1}}$ and $\beta_{M_{2}}$.

The metabolite equations detail the rates of change over time for the
intermediate substances produced and consumed by these three species, where
the concentrations are tracked for acetate in equation (2a), CO2 and H2 in
equation (2b), and polysaccharides in equation (2c). These concentrations
depend on the rate at which these metabolites are produced by the
microorganisms or enter into the system, the rate at which they are flushed
out of the system, and the rate at which they are consumed by surrounding
microorganisms. All five metabolite equations are constructed in the general
format:

$\Delta_{j}=P_{ij}-F_{j}-C_{ij},$

where $\Delta_{j}$ is the rate of change of metabolite $j$’s concentration,
$P_{ij}$ is the rate at which metabolite $j$ is produced by microorganism $i$
or enters the system, $F_{j}$ is the rate at which metabolite $j$ is flushed
out of the chemostat, and $C_{ij}$ is the rate at which the metabolite $j$ is
consumed by microorganism $i$. In each of the metabolite equations, the
parameter $\mu_{x,Y}$ is the rate at which substance $x$ is utilized by
species $Y$.

Polysaccharides enter the human gut through diet, so we accounted for their
addition to the gut through a sinusoidal function,
$\big{(}\textrm{cos}(t)+1\big{)}^{3}$. This function attempts to account for
the duration of time in between meals through the period of the curve. In
addition, this function is defined to be a smooth curve to illustrate the
gradual breakdown of food and release of nutrients in the gut.

### Parameter Estimation Techniques

In our model, many of the parameter values are unknown or cannot be determined
experimentally. In order to estimate these parameters mathematically, we
searched for data tracking the biomass changes of B. thetaiotaomicron, E.
rectale, and M smithii in order fit our model to this data with the goal of
estimating these parameters. This system, however, does not account for all
possible parameters that could potentially affect the fluctuations in biomass
or substrate concentration.

From the literature, we were able to obtain longitudinal data on an experiment
involving B. thetaiotaomicron and its substrates, presented in Adamberg et al.
2014[15]. However, similar longitudinal data on all three species B.
thetaiotaomicron, E. rectale, and M. smithii is, to our knowledge, not openly
available in the published literature. We were able to extract a single set of
data points for the full system including all three species and their relevant
substrates from Shoaie et al[11].

Due to limitations in availability of data for our 3-species model for B.
thetaiotaomicron, M. smithii, and E. rectale, we are not able to estimate the
remaining model parameters to the same degree of confidence as we were with a
reduced system (Appendix B). In our literature review, however, we were able
to obtain a single set of endpoint data values for an experiment found in
Shoaie et al.[11] conducted on all three species. Given the lack of available
longitudinal data for this biological system, we assume that this data, shown
in Table 3, are the center values of the oscillations for each substrate or
biomass quantity. One complication that arises from using this data is that
the total biomass for all three species was experimentally measured as a
single quantity, which is another factor that further contributes to our
uncertainty in our full model parameter estimates. We fit the model parameters
to produce numerical solutions with average biomass concentrations that
roughly sum to the biomass quantity given in Table 3. The parameter estimates
obtained in Appendix B on a reduced model, including some manual adjustments,
allowed us to obtain rough estimates of the remaining parameters in our full
model.

Polysaccharides ($\mu$M) | 32.06
---|---
H2 ($\mu$M) | 0
CO2 ($\mu$M) | 7.96
Acetate ($\mu$M) | 9.71
Biomass (gDW) | 0.001412
Table 3: Table of experimental data presented in Shoaie et al. 2013[11] of the microorganisms B. thetaiotaomicron, E. rectale, and M. smithii and their substrates CO2, H2, acetate, and polysaccharides. The biomass measurement is a combination of the biomass of the three microorganisms. | 0 hrs | 24 hrs | 72 hrs
---|---|---|---
Polysaccharides (mM) | 16.12 | 11.93 | 2.95
H2 (mL) | 0 | 0.0045 | 0.0051
CO2 (mL) | 0 | 0.0588 | 0.154
Acetate (mM) | 0 | 3.296 | 7.46
Biomass (gDW) | 2.34$\times$10-5 | 1.306$\times$10-4 | 2.897$\times$10-4
Table 4: Table for experimental results with data for polysaccharides, H2,
CO2, acetate, and the biomass of B. thetaiotaomicron in a medium initially
containing with 20 amino acids across distinct time points. Only the bacterial
species B. thetaiotaomicron was present in the medium. The data contained in
this table is originally from Adamberg et al. 2014 [15].

Through fitting the model parameters to the data shown in Table 3, we obtained
rough estimates of our model’s parameters, which are shown in Table 5. We
treated the parameter estimates obtained from the reduced model as initial
parameter values that could be adjusted as needed, rather than a fixed
quantity. This choice was made because the experiments conducted to produce
the data given in Tables 3 and 4 were performed under different experimental
conditions, such as the temperature of environment and the medium the culture
grew in, which are factors that could potentially impact the estimates of
these parameters.

For the remaining parameters in the model, we assigned values that were of the
same magnitude as similar parameters that we estimated in the reduced model.
After assessing the fit of the numerical solutions to the data using our
initial set of parameter estimates, these estimates were iteratively adjusted
in order to reduce the error between the center value of the numerical
solution’s oscillations and the data value given in Table 3. The resulting
parameter estimates based on implementing this approach are given in Table 5.

$\beta_{a}$ | 1.0$\times 10^{7}$ | $\beta_{h_{2}}$ | 33,000 | $\gamma_{B}$ | 10 | $\mu_{hM}$ | 40
---|---|---|---|---|---|---|---
$\beta_{B}$ | 1.2 | $\beta_{M_{1}}$ | 0.8 | $\gamma_{h}$ | 150 | $\mu_{pB}$ | 200,000
$\beta_{E_{1}}$ | 0.8 | $\beta_{M_{2}}$ | 0.5 | $\gamma_{p}$ | 400 | $\mu_{pE}$ | 5,000
$\beta_{E_{2}}$ | 0.6 | $\beta_{p}$ | 10,000 | $\mu_{aE}$ | 25,000 | $q$ | 0.05
$\beta_{h_{1}}$ | 150 | $\gamma_{a}$ | 200 | $\mu_{aM}$ | 50,000 |  |
Table 5: Table of fitted parameter values for equations (1) and (2) based on
the experimental data in Table 4.

### Sensitivity Analysis Methods

We utilized global sensitivity analysis in order to identify the driving
parameters of our model. In this framework, which utilizes Monte-Carlo (MC)
methods, each parameter is assigned a probability density function (pdf) based
on a priori information known about the parameter. Samples are drawn from the
probability density functions (pdfs) to evaluate the overall model output.

Because we are interested in detecting significant nonlinear interactions
among the model parameters, implementations of global sensitivity analysis
methods for linear relationships like the partial rank correlation coefficient
(PRCC), the Pearson correlation coefficient (CC), and standardized regression
coefficients (SRC) would not be useful in our case. Consequently, we employed
the Sobol’ method, a variance-based decomposition method[26].

Sobol’ indices are a MC variance-based approach to calculating all first-order
and total-effects indices in a model with $k$ parameters. These sensitivity
indices are computed based on model evaluations for $N$ simulations. In order
to improve the sensitivity estimates, the values tested for each parameter in
the model are drawn from quasi-random number generators. Though the use of
quasi-random numbers from a given distribution is not necessary, we
incorporate this approach into our computations.

We wrote a function in Julia v1.1.1 that computes the first-order and total-
order Sobol’ indices for a given ODE-based model for each parameter-variable
combination, which is included in Appendix A.

| $B$ | E | $M$ | $a$ | $h$ | $p$
---|---|---|---|---|---|---
$\beta_{a}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\beta_{B}$ | $0.98$ | $0.02$ | $0.01$ | $0.02$ | $0.37$ | $0.13$
$\beta_{E_{1}}$ | $0.00$ | $0.00$ | $0.00$ | $0.04$ | $0.00$ | $0.00$
$\beta_{E_{2}}$ | $0.00$ | $0.82$ | $0.01$ | $0.02$ | $0.10$ | $0.03$
$\beta_{h_{1}}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\beta_{h_{2}}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\beta_{M_{1}}$ | $0.00$ | $0.00$ | $0.04$ | $0.04$ | $0.06$ | $0.00$
$\beta_{M_{2}}$ | $0.00$ | $0.00$ | $0.65$ | $0.02$ | $0.17$ | $0.00$
$\beta_{p}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\gamma_{a}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\gamma_{B}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\gamma_{h}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\gamma_{p}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.01$ | $0.00$
$\mu_{aE}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\mu_{aM}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\mu_{hM}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\mu_{pB}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\mu_{pE}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$q$ | $0.00$ | $0.01$ | $0.02$ | $0.07$ | $0.02$ | $0.09$
Total | $0.99$ | $0.85$ | $0.74$ | $0.23$ | $0.75$ | $0.24$

Table 6: Table of model parameters’ estimated first-order Sobol’ indices for
each output, $B$, $E$, $M$, $a$, $h$, and $p$, using $N=2^{18}$ simulations.

| $B$ | E | $M$ | $a$ | $h$ | $p$
---|---|---|---|---|---|---
$\beta_{a}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\beta_{B}$ | $0.99$ | $0.13$ | $0.11$ | $0.33$ | $0.56$ | $0.80$
$\beta_{E_{1}}$ | $0.00$ | $0.01$ | $0.00$ | $0.56$ | $0.01$ | $0.04$
$\beta_{E_{2}}$ | $0.01$ | $0.97$ | $0.04$ | $0.22$ | $0.18$ | $0.49$
$\beta_{h_{1}}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\beta_{h_{2}}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\beta_{M_{1}}$ | $0.00$ | $0.00$ | $0.05$ | $0.58$ | $0.10$ | $0.01$
$\beta_{M_{2}}$ | $0.00$ | $0.00$ | $0.88$ | $0.19$ | $0.25$ | $0.00$
$\beta_{p}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\gamma_{a}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\gamma_{B}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\gamma_{h}$ | $0.00$ | $0.00$ | $0.01$ | $0.00$ | $0.00$ | $0.00$
$\gamma_{p}$ | $0.01$ | $0.01$ | $0.00$ | $0.00$ | $0.01$ | $0.00$
$\mu_{aE}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\mu_{aM}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\mu_{hM}$ | $0.00$ | $0.00$ | $0.04$ | $0.01$ | $0.00$ | $0.00$
$\mu_{pB}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$\mu_{pE}$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$q$ | $0.01$ | $0.05$ | $0.11$ | $0.62$ | $0.17$ | $0.72$
Total | $1.02$ | $1.17$ | $1.25$ | $2.48$ | $1.29$ | $2.09$

Table 7: Table of model parameters’ estimated total-order Sobol’ indices for
each output, $B$, $E$, $M$, $a$, $h$, and $p$, using $N=2^{18}$ simulations.

To compute the first- and total-order Sobol’ indices for our reduced model
parameters, we specified prior distributions for each parameter from which to
draw samples. Parameter samples are initially drawn from a quasi-random number
generator on the interval [0,1]. Using the inverse transform sampling method,
these samples are transformed into a simulated parameter value based on its
specified prior distribution for the sensitivity calculation of the Sobol’
indices[27].

To specify distributions for each parameter in the full model, we first set
all but one parameter to a Gamma($k,\theta$) distribution since we expect all
of these parameters to be strictly positive. The parameter $q$ was specified
to follow a Beta$(5,95)$ distribution. Though the support of this distribution
restricts the samples to the interval $[0,1]$, we deemed this option to be
more appropriate due to the Gamma(0.025, 2.0) distribution’s heavy right-
skewness and tendency to sample infinitesimally small, unlikely values.

We implemented a similar approach in specifying each $k$ and $\theta$ as in
the reduced model; the mean of each parameter’s distribution was set to be
equal to our estimated value in Table 5. The variance was specified to be
equal to two times the estimated parameter value. We chose this value for the
variance in order to cover a wide range of sample values for our sensitivity
analysis due to a large uncertainty in our parameter estimates. Table 1 shows
the values of $k$, $\theta$, and the mean and variance of each parameter
distribution.

## Data availability

The data generated in this modeling and simulation effort is included in the
appendices of this manuscript, as well as the Julia code used to generate that
data.

## References

  * [1] David, L. A. _et al._ Diet rapidly and reproducibly alters the human gut microbiome.  _Nature_ 505, 559–563, DOI: 10.1038/nature12820.Diet (2014).
  * [2] Ji, B. & Nielsen, J.  From next-generation sequencing to systematic modeling of the gut microbiome.  _Frontiers in Genetics_ 6, 1–9, DOI: 10.3389/fgene.2015.00219 (2015).
  * [3] Mangalam, A. K., Yadav, M. & Yadav, R.  The emerging world of microbiome in autoimmune disorders: Opportunities and challenges.  _Indian journal of rheumatology_ 16, 57—72, DOI: 10.4103/injr.injr_210_20 (2021).
  * [4] Freedman, S. N., Shahi, S. K. & Mangalam, A. K.  The “gut feeling”: Breaking down the role of gut microbiome in multiple sclerosis.  _Neurotherapeutics_ 1, 109—125, DOI: 10.1007/s13311-017-0588-x (2018).
  * [5] Cryan, J. F., O’Riordan, K. J., Sandhu, K., Peterson, V. & Dinan, T. G.  The gut microbiome in neurological disorders.  _Lancet Neurol._ 2, 179—194, DOI: 10.1016/S1474-4422(19)30356-4 (2020).
  * [6] Singh, R. K. _et al._ Influence of diet on the gut microbiome and implications for human health.  _J Transl Med._ 1, 73, DOI: 10.1186/s12967-017-1175-y (2017).
  * [7] Cady, N., Peterson, S. R., Freedman, S. N. & Mangalam, A. K.  Beyond metabolism: The complex interplay between dietary phytoestrogens, gut bacteria, and cells of nervous and immune systems.  _Frontiers in Neurology_ 11, 150, DOI: 10.3389/fneur.2020.00150 (2020).
  * [8] Kumar, M., Ji, B., Zengler, K. & Nielsen, J.  Modelling approaches for studying the microbiome.  _Nature Microbiology_ 4, 1253–1267, DOI: 10.1038/s41564-019-0491-9 (2019).
  * [9] Oreskes, N., Shrader-Frechette, K. & Belitz, K.  Verification, validation, and confirmation of numerical models in the earth sciences.  _Science_ 263, 641–646, DOI: 10.1126/science.263.5147.641 (1994).
  * [10] Smith, H. L. & Waltman, P.  _The Theory of The Chemostat: Dynamics of Microbial Competition_ (Cambridge University Press, 2008).
  * [11] Shoaie, S. _et al._ Understanding the interactions between bacteria in the human gut through metabolic modeling.  _Scientific Reports_ 3, 1–10, DOI: 10.1038/srep02532 (2013).
  * [12] Kim, M. H. & Kim, H.  The roles of glutamine in the intestine and its implication in intestinal diseases.  _International Journal of Molecular Sciences_ 18, DOI: 10.3390/ijms18051051 (2017).
  * [13] Scaldaferri, F. _et al._ Intestinal gas production and gastrointestinal symptoms: from pathogenesis to clinical implication.  _European review for medical and pharmacological sciences_ 17, 2–10 (2013).
  * [14] Xu, J. _et al._ A genomic view of the human-Bacteroides thetaiotaomicron symbiosis.  _Science_ 299, 2074–2076, DOI: 10.1210/jcem-10-10-1361 (2003).
  * [15] Adamberg, S. _et al._ Degradation of fructans and production of propionic acid by Bacteroides thetaiotaomicron are enhanced by the shortage of amino acids.  _Frontiers in Nutrition_ 1, 1–10, DOI: 10.3389/fnut.2014.00021 (2014).
  * [16] Martens, E. C. _et al._ Recognition and degradation of plant cell wall polysaccharides by two human gut symbionts.  _PLoS Biology_ 9, DOI: 10.1371/journal.pbio.1001221 (2011).
  * [17] Glass, T. L. & Hylemon, P. B.  Characterization of a pyridine nucleotide-nonspecific glutamate dehydrogenase from bacteroides thetaiotaomicron.  _Journal of Bacteriology_ 141, 1320–1330 (1980).
  * [18] Mahowald, M. A. _et al._ Characterizing a model human gut microbiota composed of members of its two dominant bacterial phyla.  _Proceedings of the National Academy of Sciences of the United States of America_ 106, 5859–5864 (2009).
  * [19] Saltelli, A. _et al._ _Global Sensitivity Analysis. The Primer_ , chap. 1-6, i–xi (John Wiley & Sons, Ltd, 2008).  https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470725184.fmatter.
  * [20] Iooss, B. & Prieur, C.  Shapley effects for sensitivity analysis with correlated inputs: Comparisons with sobol’ indices, numerical estimation and applications.  _International Journal for Uncertainty Quantification_ 9, 493–514, DOI: 10.1615/Int.J.UncertaintyQuantification.2019028372 (2019).  1707.01334.
  * [21] Rabitz, H.  Sixth International Conference on Sensitivity Analysis of Model Output Global Sensitivity Analysis for Systems with Independent and / or Correlated Inputs.  _Soins Aides-soignantes_ 2, 7587–7589, DOI: 10.1016/j.sbspro.2010.05.131 (2010).
  * [22] Jangi, S. _et al._ Alterations of the human gut microbiome in multiple sclerosis.  _Nature Communications_ 7, DOI: 10.1038/ncomms12015 (2016).
  * [23] Monod, J.  The growth of bacterial cultures.  _Ann. Rev. Microbiol._ 3, 371–394 (1949).
  * [24] Stefan, M. I. & Le Novère, N.  Cooperative Binding.  _PLoS Computational Biology_ 9, e1003106, DOI: 10.1371/journal.pcbi.1003106 (2013).
  * [25] Voit, E. O.  _Computational Analysis of Biochemical Systems: A Practical Guide for Biochemists and Molecular Biologists_ (Cambridge University Press, 2000).
  * [26] Marino, S., Hogue, I. B., Ray, C. J. & Kirschner, D. E.  A methodology for performing global uncertainty and sensitivity analysis in systems biology.  _Microbiology_ 254, 178–196, DOI: 10.1016/j.jtbi.2008.04.011.A (2009).
  * [27] Devroye, L.  General principles in random variate generation.  In _Non-Uniform Random Variate Generation_ , chap. 2, 27–82 (Springer, 1986).

## Author contributions statement

M.A. translated the biological system to a schematic, wrote the Julia code,
and conducted the numerical simulations and sensitivity analyses. M.A and B.A.
created the mathematical models. M.A. and A.M. related the models and
underlying biology to each other. B.A. conceived of the project. All authors
contributed to the writing of the manuscript and reviewed the manuscript.

## Additional information

Competing interests

The authors declare no competing interests.

##  Appendix A Appendix A: Julia Code for First- and Total-Order Sobol’
Indices

⬇

"""

sobol(f,␣N,␣u0,␣tspan,␣t,␣k,␣$\theta$,␣index␣=␣"first-order")

Creates␣a␣matrix␣of␣Sobol’␣sensitivities␣where␣each␣row␣contains

the␣sensitivity␣indices␣for␣a␣parameter␣across␣different␣model

outputs.

␣f␣-␣system␣of␣ODEs

␣N␣-␣number␣of␣simulations␣(can␣vary␣from␣hundreds␣to␣thousands)

␣u0␣-␣vector␣of␣initial␣conditions␣for␣ODE

␣tspan␣-␣time␣span

␣k␣-␣vector␣of␣scale␣parameters␣from␣Gamma($k_{i}$,$\theta_{i}$)␣for␣all

␣␣␣␣␣model␣parameters

␣$\theta$␣-␣vector␣of␣shape␣parameters␣from␣Gamma($k_{i}$,$\theta_{i}$)␣for␣all

␣␣␣␣␣model␣parameters

␣index␣-␣default␣is␣the␣first␣order␣index,␣but␣can␣take␣on␣the

␣␣␣␣␣␣␣␣␣other␣value␣"total-order"

␣t␣-␣specified␣time␣in␣tspan␣to␣evaluate␣each␣model␣variable

"""

using Sobol

using Distributions

function sobol(f, N, u0, tspan, k, $\theta$, index = "first-order",

t = 5)

K = length(k)

vars = length(u0)

# quasi-randomly selected points

S = zeros(N,2*K)

s = SobolSeq(2*K)

$\theta$ = [$\theta$;$\theta$]

k = [k;k]

S = zeros(N,2*K)

for sim in 1:N

S[sim,:] = next!(s)”

⬇

for parameter in 1:2*K

S[sim,parameter] = quantile(Gamma(k[parameter],

$\theta$[parameter]), S[sim,parameter])

end

end

A = S[:, 1:K] # accounts for half of the random sample

B = S[:, (K+1):2*K] # accounts for second half of the

# random sample

sums = zeros(4,vars)

ya = solveODE(f, u0, vars, tspan, t, A, N)

yb = solveODE(f, u0, vars, tspan, t, B, N)

sumC = zeros(K,vars) # (parameter, variable)

for p in 1:K # iterates through changing one column of C at

# a time (parameters)

C = Ci(p,A,B)

yc = solveODE(f, u0, vars, tspan, t, C, N)

for v in 1:vars

if index == "first-order"

sumC[p,v] = yc[:,v]’*ya[:,v]/N’

⬇

else

sumC[p,v] = yc[:,v]’*yb[:,v]/N’

⬇

end

end

end

sens = zeros(K,vars) # (parameter, variable)

⬇

for v in 1:vars

sums[1,v] = ya[:,v]’*ya[:,v]/N’

⬇

# sum(ya^2)/N

⬇

sums[2,v] = (sum(ya[:,v])/N)^2 # [sum(ya)/N]^2

sums[3,v] = sums[1, v] - sums[2,v]

# denominator for sensitivity calculation

sums[4,v] = sum(yb[:,v])*sum(ya[:,v])/N.^2

# sum(ya)*sum(yb)/N^2

for p in 1:K

if index == "first-order"

sens[p,v] = (sumC[p,v] - sums[4,v])/sums[3,v]

else

sens[p,v] = 1 - (sumC[p,v] - sums[2,v])/sums[3,v]

end

end

end

return sens # matrix of sensitivity indices

end

⬇

"""

Ci(i,␣A,␣B)

Creates␣a␣new␣matrix␣C␣by␣replacing␣the␣ith␣column␣of␣matrix␣B

with␣the␣ith␣column␣of␣matrix␣A.

␣i␣-␣column␣number␣to␣be␣changed␣in␣matrix␣B

␣A␣-␣matrix␣containing␣half␣of␣the␣pseudo-random␣samples

␣B␣-␣matrix␣containing␣the␣other␣half␣of␣the␣pseudo-random

␣␣␣␣␣samples

"""

function Ci(i, A, B)

C = zeros(size(B))

for j in 1:size(A)[2]

if j == i

C[:,i] = A[:,i]

else

C[:,j] = B[:,j]

end

end

return C

end

⬇

"""

solveODE(f,␣u0,␣vars,␣tspan,␣t,␣M,␣N)

Solves␣the␣system␣of␣ODEs␣f␣and␣creates␣a␣matrix␣of␣model

evaluations␣for␣each␣simulation␣and␣output␣variable.␣Each␣row

corresponds␣to␣model␣solutions␣at␣a␣given␣time␣t␣for␣one

simulation.

␣f␣-␣system␣of␣ODEs

␣u0␣-␣vector␣of␣initial␣conditions␣for␣the␣ODEs

␣vars␣-␣number␣of␣output␣variables

␣tspan␣-␣time␣span

␣t␣-␣specified␣time␣in␣tspan␣to␣evaluate␣each␣model␣variable

␣M␣-␣matrix␣with␣rows␣as␣parameter␣values␣and␣columns␣as

␣␣␣␣␣simulations

␣N␣-␣number␣of␣simulations␣(can␣vary␣from␣hundreds␣to␣thousands)

"""

function solveODE(f, u0, vars, tspan, t, M, N)

y = zeros(N,vars) # (simulation, variable)

for sim in 1:N

p = M[sim,:]

prob = ODEProblem(f,u0,tspan,p)

sol = solve(prob)

for variable in 1:vars

y[sim,variable] = sol.u[t][variable]

end

end

return y # matrix of model outputs

end

##  Appendix B Appendix B: Reduced System

Table 1 presents data found in Adamberg et al. 2014[15]. The experiment
conducted in this paper focused on B. thetaiotaomicron with the substrates H2,
CO2, acetate, and polysaccharides. We utilized this data for parameter
determination and estimation by condensing our model to only include B.
thetaiotaomicron and its substrates.

| 0 hrs | 24 hrs | 72 hrs
---|---|---|---
Polysaccharides (mM) | 16.12 | 11.93 | 2.95
H2 (mL) | 0 | 0.0045 | 0.0051
CO2 (mL) | 0 | 0.0588 | 0.154
Acetate (mM) | 0 | 3.296 | 7.46
Biomass (gDW) | 2.34$\times$10-5 | 1.306$\times$10-4 | 2.897$\times$10-4
Table 8: Table for experimental results with data for polysaccharides, H2,
CO2, acetate, and the biomass of B. thetaiotaomicron in a medium initially
containing with 20 amino acids across distinct time points. Only the bacterial
species B. thetaiotaomicron was present in the medium. The data contained in
this table is originally from Adamberg et al. 2014 [15].

The subsystem isolating B. thetaiotaomicron’s metabolism is illustrated in
Figure 6. Based on this graphical representation, we extracted the reduced
system from our full model in equations (1) and (2), producing the reduced
system of ODEs in equations (5). Terms related to chemostatic flow were
removed because the experimental data in Table 8 was collected in a closed
system, rather than an open system like a chemostat. Lastly, we removed any
terms that account for an input of polysaccharides because only a specified
initial amount was placed in the experimental device for B. thetaiotaomicron
to consume over time.

Figure 6: Graphical representation of a subset of B. thetaiotaomicron’s
metabolism, including the utilization of polysaccharides and the production of
acetate, CO2, and H2. This diagram is based on information presented in
Adamberg et al. 2014[15].

$\displaystyle\frac{dB}{dt}$ $\displaystyle=\beta_{B}\Psi_{\gamma_{p}}(p)B$
(5a) $\displaystyle\frac{da}{dt}$
$\displaystyle=\beta_{a}\Psi_{\gamma_{p}}(p)B$ (5b)
$\displaystyle\frac{dh}{dt}$
$\displaystyle=\beta_{h_{2}}\Psi_{\gamma_{p}}(p)B$ (5c)
$\displaystyle\frac{dp}{dt}$ $\displaystyle=-\mu_{p,B}Bp$ (5d)

By fitting our reduced model to the collected data in Table 8, we estimated
the model parameters using the Julia packages `DiffEqParamEstim` and `Optim`
with the goal of utilizing these results to inform our full model parameter
estimation. We ultimately chose the simulated annealing method to apply to our
optimization problem.

In the first step of finding a roughly optimal solution for our reduced
model’s parameters, we first chose an arbitrary vector of initial values for
these parameters. Once we chose these arbitrary initial values, we used the
`DiffEqParamEstim` and `Optim` packages in order to obtain the optimal results
given by the simulated annealing method. The results produced from this
constrained optimization given our initial guess of parameter values allowed
us to achieve a seemingly local optimum for the system. We utilized this
output by adjusting the output parameter values ourselves in order to further
reduce the residuals. Using this new set of parameters, we reran the
optimization once more in the hopes of producing a better fit to the data.

Figures 7, 8, and 9 show the solutions to our system of ODEs using our final
set of roughly optimized parameters in Table 9 based on the data in Table 8.

Figure 7: Plot of the solutions to equation (5a) and the observed data for B. thetaiotaomicron over the time interval 0 to 72 hours. Figure 8: Plot of the solutions to equations (5b) and (5d), as well as the observed data for acetate and polysaccharides over the time interval 0 to 72 hours. Figure 9: Plot of the solutions to equation (5c) and the observed data for CO2 and H2 over the time interval 0 to 72 hours. $\beta_{a}$ | $\beta_{B}$ | $\beta_{h_{2}}$ | $\gamma_{p}$ | $\mu_{pB}$
---|---|---|---|---
2000 | 0.08 | 50 | 10 | 200
Table 9: Table of fitted parameter values for equations (5) based on the
experimental data in Table 8.

#### First- and Total-Order Sobol’ Indices for the Reduced Model Parameters

In order to specify the distribution of our model’s parameters, we first
considered the domain of each parameter. All parameters in our reduced model
can reasonably take on positive values, so we restricted our consideration to
parameter distributions defined on the positive number line. Due to its
positive, continuous domain, we chose the Gamma($k$, $\theta$) distribution as
the general form of our specified parameter distributions.Using the results of
our parameter estimation in Table 9, we defined the mean of the Gamma
distributions to be equal to the estimates we obtained from our parameter
estimation. Because we obtained one set of parameter estimates for our reduced
model, we are not able to obtain an informed variance on our parameter
distributions. In order to reflect our minimal knowledge of the true parameter
values, we define the variance of these Gamma distributions to be large in
order to test a wide range of parameter values in our sensitivity
calculations, which are specified in Table 10.

| $k$ | $\theta$ | Mean | Variance
---|---|---|---|---
$\beta_{a}$ | 2000 | 1 | 2000 | 2000
$\beta_{B}$ | 0.01 | 8 | 0.08 | 0.64
$\beta_{h}$ | 50 | 1 | 50 | 50
$\gamma_{p}$ | 10 | 1 | 10 | 10
$\mu_{pB}$ | 200 | 1 | 200 | 200
Table 10: Table of defined parameter distributions for the parameters in our reduced model, including the mean and variance of the distribution. | $B$ | $a$ | $h$ | $p$
---|---|---|---|---
$\beta_{a}$ | 0.000 | 0.003 | 0.002 | 0.003
$\beta_{B}$ | 0.995 | 0.838 | 0.572 | 0.891
$\beta_{h}$ | 0.000 | 0.000 | 0.304 | 0.000
$\gamma_{p}$ | 0.000 | 0.120 | 0.080 | 0.064
$\mu_{pB}$ | 0.000 | 0.025 | 0.017 | 0.034
Total | 0.995 | 0.982 | 0.974 | 0.992
Table 11: Table of first-order Sobol’ indices for our reduced model in equations (5) using $N=2^{18}$ simulations. | $B$ | $a$ | $h$ | $p$
---|---|---|---|---
$\beta_{a}$ | 0.003 | 0.003 | 0.002 | 0.002
$\beta_{B}$ | 1.000 | 0.856 | 0.596 | 0.898
$\beta_{h}$ | 0.002 | 0.001 | 0.318 | 0.000
$\gamma_{p}$ | 0.003 | 0.133 | 0.092 | 0.072
$\mu_{pB}$ | 0.002 | 0.030 | 0.021 | 0.035
Total | 1.010 | 1.012 | 1.028 | 1.007
Table 12: Table of total-order Sobol’ indices for our reduced model in
equations (5) using $N=2^{18}$ simulations.

The first- and total-order results of the sensitivity analysis of the
parameters in our reduced model are presented in Tables 11 and 12. Based on
the small changes observed from the first-order sensitivities compared to the
total-order sensitivities, we can conclude that higher order parameter
interactions do not significantly contribute to the variance of the model
outputs, so these higher order interactions can safely be ignored. From the
sensitivity results, it can be concluded that the parameter $\beta_{B}$ is
influential on the variance of the model outputs due to $\beta_{B}$’s high
sensitivity indices for all outputs. An unexpected result based on these
indices is that the variable $p$ is more sensitive to changes in the parameter
$\gamma_{p}$, the Michaelis-Menten constant in the Monod form, than
$\mu_{pB}$, the rate B. thetaiotaomicron consumes polysaccharides. Though
$\gamma_{p}$ explains more variance in $p$ than $\mu_{pB}$, the first- and
total-order effects for this parameter are relatively low compared to the main
effects of parameter $\beta_{B}$, which appears to be the driving parameter of
this subsystem. Overall, the results obtained in Tables 11 and 12 are
generally unsurprising and intuitive. With these parameter estimates and
Sobol’ indices in mind, we then utilized this information in our full model to
conduct a sensitivity analysis of all nineteen parameters.