of 11
Khazaei
et al
.,
Sci. Adv.
2020;
6
: eaba0353 12 August 2020
SCIENCE ADVANCES
|
RESEARCH ARTICLE
1 of 10
MICROBIOLOGY
Metabolic multistability and
hysteresis in
a
model
aerobe-anaerobe microbiome community
Tahmineh
Khazaei
1
, Rory L.
Williams
1
, Said R.
Bogatyrev
1
, John C.
Doyle
1,2
,
Christopher S.
Henry
3
, Rustem F.
Ismagilov
1,4
*
Major changes in the microbiome are associated with health and disease. Some microbiome states persist despite
seemingly unfavorable conditions, such as the proliferation of aerobe-anaerobe communities in oxygen-exposed
environments in wound infections or small intestinal bacterial overgrowth. Mechanisms underlying transitions
into and persistence of these states remain unclear. Using two microbial taxa relevant to the human microbiome,
we combine genome-scale mathematical modeling, bioreactor experiments, transcriptomics, and dynamical sys-
tems theory to show that multistability and hysteresis (MSH) is a mechanism describing the shift from an aerobe-
dominated state to a resilient,
paradoxically persistent aerobe-anaerobe state. We examine the impact of changing
oxygen and nutrient regimes and identify changes in metabolism and gene expression that lead to MSH and as-
sociated multi-stable states. In such systems, conceptual causation-correlation connections break and MSH must
be used for analysis. Using MSH to analyze microbiome dynamics will improve our conceptual understanding of
stability of microbiome states and transitions between states.
INTRODUCTION
Recent evidence shows that changes in the species composition and
abundance of the human microbiome can be associated with health
and disease (
1
3
). Understanding the mechanisms that cause com-
positional shifts in healthy microbiomes, which otherwise can be stable,
is challenging because of the inherent complexity of these ecosystems.
A perplexing feature of some of these disturbed ecosystems is the
persistence of a new microbiome state, even in seemingly unfavorable
conditions. For example, in small intestinal bacterial overgrowth (SIBO),
strict anaerobes that are typically found only in the colon become
prominent in the small intestine and, paradoxically, persist in this
environment exposed to oxygen flux from the tissue (
4
,
5
). Similar-
ly, in periodontal diseases (
6
) and in wound infections, anaerobes
proliferate in oxygen-exposed environments.
One potential mechanism to explain microbiome shifts and their
persistence is multistability (
7
13
), the concept that several steady
states can exist for an identical set of system parameters (Fig. 1).
Multistable systems have been described in the context of ecosys
-
tems (
14
17
) and gene regulatory networks (
18
20
). Now, with the
expanding characterization of the microbiome, there are signs that
multistability may also exist in these communities (
21
27
). For ex-
ample, compositional changes in gut microbiota are implicated in
inflammatory bowel disease (
28
) and obesity (
29
). Bimodal species
abundance (i.e., when a microbial species is present at either high or
low levels) has been interpreted as multistability (
30
); however, as
discussed by Gonze
et al.
, bimodality is insufficient to prove mul-
tistability (
7
,
31
). Some multistable systems can additionally exhibit
hysteresis, where in response to a perturbation, a system gets “stuck” in
a new steady state and the former state cannot be regained by simply
reversing the perturbation (
7
). The presence of hysteresis could be
hypothesized from studies of the microbiome (
32
). For example, antibiotic
exposures can change the microbiome composition and have lasting
effects even after removal of the antibiotic (
33
,
34
). However, it has not
been rigorously tested whether multistability and hysteresis (MSH) can
arise in a microbiome-relevant community and by what mechanism.
Here, we investigate MSH in a minimally “complex” two-species
system to represent the paradoxical aerobe-anaerobe microbiome com-
munities that persist in oxygen-exposed environments. We used two
organisms prevalent in SIBO (
35
): the anaerobe
Bacteroides thetaio-
taomicron
(
Bt
) that breaks down complex carbohydrates (e.g., dex-
tran) into simple sugars and short-chain fatty acids (
36
) and the
facultative anaerobe (hereafter referred to as an aerobe)
Klebsiella
pneumoniae
(
Kp
) capable of consuming oxygen, simple sugars, and
short-chain fatty acids and performing anaerobic respiration in the
absence of oxygen (
37
).
Mathematical simulations (Fig. 2A) and a 35-day (832 hours)
experiment (fig. S1) in a continuously stirred tank reactor (CSTR)
revealed that MSH occurs in our two-species system with two distinct
steady states that can exist under identical environmental conditions:
(i)
Kp
-only state and (ii)
Kp-Bt
state, where the anaerobe (
Bt
) prolifer-
ates in the presence of continuous oxygen input. Using genome-scale
mathematical models, which capture the full metabolic capacity of
each species (1491 equations for
Bt
and 2262 equations for
Kp
), and
RNA sequencing data collected from the CSTR experiment, we find
that MSH extends to the level of metabolism, where genes are dif-
ferentially expressed in the two distinct states. We identified key
metabolic pathways (short-chain fatty acid and oligosaccharide me-
tabolism) involved in metabolic coupling between the two species
leading to MSH.
RESULTS
Computational results
Mathematical simulations of a CSTR revealed that MSH arises from
the interplay between environmental perturbations and interspecies
metabolic interactions. We used the dynamic multispecies metabolic
modeling (DMMM) framework (
38
) to model a community of
Kp
1
Division of Biology and Biological Engineering, California Institute of Technology,
Pasadena, CA, USA.
2
Division of Engineering and Applied Science, California Insti-
tute of Technology, Pasadena, CA, USA.
3
Data Science and Learning Division,
Argonne National Laboratory, Lemont, IL, USA.
4
Division of Chemistry and Chemical
Engineering, California Institute of Technology, Pasadena, CA, USA.
*Corresponding author. Email: rustem.admin@caltech.edu
Copyright © 2020
The
Authors, some
rights reserved;
exclusive licensee
American Association
for the Advancement
of Science. No claim to
original U.S.
Government
Works. Distributed
under a Creative
Commons Attribution
License 4.0 (CC BY).
on August 13, 2020
http://advances.sciencemag.org/
Downloaded from
Khazaei
et al
.,
Sci. Adv.
2020;
6
: eaba0353 12 August 2020
SCIENCE ADVANCES
|
RESEARCH ARTICLE
2 of 10
and
Bt
in a CSTR (Fig. 2A) with continuous input flows of dextran
minimal media and varying input glucose or varying input oxygen
levels (depending on the simulation). The outflow rate is equal to the
inflow rate to maintain a constant reactor volume, and the resident
time in the reactor is 5 hours. The DMMM framework uses dynamic
flux balance analysis (dFBA) (
39
), which allows us to capture tem-
poral changes in intracellular flux rates (using the genome-scale
metabolic model for each species), extracellular metabolite concen-
trations, and species concentrations.
To computationally test whether a change in the balance of oxygen
and carbon fluxes could lead to a change in the state of the aerobe-
anaerobe community, we altered glucose input concentrations (Fig. 3A),
while keeping constant all other system parameters, including con-
tinuous oxygen input and continuous dextran input. The model
predicted that for glucose concentrations of 0.25 to 3 mM in the
input feed (at a constant flow rate of 0.7 ml/min for all conditions), the
output state consisted solely of
Kp
, which we refer to as the
Kp
-only
state (Fig. 2B). Stoichiometrically, at these glucose concentrations,
oxygen was not completely consumed; thus, the environment was
unfavorable for
Bt
growth. However, when we increased glucose in-
put concentration to 3.25 mM, we observed a shift to a new steady
state (Fig. 3A). At this “tipping point,” the environment became suffi-
ciently anaerobic to support the growth of
Bt
. We refer to this second
distinct steady state as the
Kp-Bt
(aerobe-anaerobe) state (Fig. 2C).
In the
Kp-Bt
state, the
Kp
population is no longer carbon limited
because of the additional carbon sources generated from the metab-
olism of dextran by
Bt
. Thus,
Kp
can now consume all of the available
oxygen to oxidize both glucose and the additional carbon sources,
resulting in anaerobic conditions. Unexpectedly, this
Kp-Bt
state
persisted when we systematically reversed the input of glucose be-
low 3.25 mM, even to 0 mM.
Thus, this system shows hysteresis and
multistability: Under identical input conditions of glucose and oxy
-
gen, the system can be in either of the two possible states. We then
identified tipping points for population shifts in response to input oxy
-
gen variations, with glucose kept constant (Fig. 3B). In addition to
glucose input, all other parameters, including continuous dextran
input, were held constant. We considered oxygen as a parameter because
in host settings, oxygen availability can be affected by respiration,
blood flow rate, immune consumption, etc. We found that we could
return the system to the
Kp
-only state by increasing oxygen levels, a
state switch that was not possible by manipulating glucose concen-
tration alone. Last, we simulated changes in both glucose and oxygen
levels and characterized the landscape of multistability and monosta
-
bility in the model microbial community (Fig. 3C). These simulation
results illustrate that even a minimal model of microbiome with code-
pendence (
40
,
41
) can demonstrate marked MSH.
Experimental results
To confirm and further explore MSH beyond mathematical predic-
tions, we performed a CSTR experiment over 35 days (832 hours).
In the CSTR (200-ml culture volume), we varied input glucose con-
centrations (while keeping all other parameters constant, including
continuous dextran and oxygen) and measured the steady-state out-
put composition of the microbial community by quantitative poly-
merase chain reaction [qPCR; and digital PCR (dPCR); fig. S6]. Oxygen
was sparged into the reactor at 3.4% of the gas feed (total gas feed,
50 ml/min) and kept constant for all conditions. For each steady-state
condition, we collected three CSTR samples separated by at least one
residence time (5 hours; fig. S1 contains the experimental workflow).
As predicted by the mathematical models, we observed both mul-
tistability and hysteresis (Fig. 4A) experimentally. At 0.25, 1, and 2 mM
glucose concentrations, the steady-state community consisted only of
Kp
;
Bt
was washed out under these conditions (Fig. 4B). To confirm
washout,
Bt
was reinoculated three separate times. The dissolved oxy-
gen measurements (Fig. 4C) confirmed that oxygen was not limiting
under the selected parameter conditions, resulting in an aerobic en-
vironment unsuitable for
Bt
growth. As in the simulations, at 5 mM
glucose, a new distinct steady state was reached where
Bt
grew in
the presence of
Kp.
Although there was continuous oxygen flux
into the reactor, the concentration of dissolved oxygen measured in
the reactor was near zero. Next, to test for hysteresis, we reduced the
glucose input back down to 2, 1, 0.25, and 0 mM and found that the
aerobe-anaerobe state persisted. The persistence of the
Kp-Bt
state
(instead of a return to the
Kp
-only state) qualitatively confirmed model
predictions of hysteresis and verified that this microbial community
is a multistable system.
The CSTR results demonstrate metabolic coupling and code-
pendence between these two bacterial species with respect to carbon
and oxygen. At sample point 8, there is no glucose input to the reactor,
yet
Kp
continued to grow, indicating that
Kp
was completely depen-
dent on
Bt
for its carbon supply. At sample point 4,
Bt
started to grow,
Seesaw state
Left state
Right state
Frictionless seesaw
Seesaw state
Distance
Left state
Right state
“Tipping point”
“Tipping point”
Hysteresis
1
1
2
3
4
2
3
4
1
8
4
5
2
7
3
6
1
2
3
4
5
6
7
8
“Tipping point”
Seesaw with friction
Fig. 1. Simplified illustration of multistability and hysteresis.
A frictionless see-
saw is a multistable system with two states (left and right) determined by the posi-
tion of the person. Conversely, a seesaw with static friction is a multistable system
with hysteresis. Within the region of hysteresis, the state of the system is not deter-
mined by the external input (position of person) but rather by the system’s history.
As the person walks through positions 1 and 2 past the midpoint to position 3, the
seesaw remains stuck in the left state; the person has to walk much further, to po-
sition 4, to switch it to the right state. As the person walks back to the left, through
positions 5 and 6 and past the midpoint to position 7, the seesaw remains stuck in
the right state and it only switches back to the left state when the person reaches
position 8. In the region of positions 2, 3, 6, and 7, the system exhibits multistability
and hysteresis (MSH). Here, the state of the microbial community (
Kp
-only state or
Kp-Bt
state) is analogous to the state of the seesaw, and the balance of oxygen or
glucose inputs are analogous to the positions of the person.
on August 13, 2020
http://advances.sciencemag.org/
Downloaded from
Khazaei
et al
.,
Sci. Adv.
2020;
6
: eaba0353 12 August 2020
SCIENCE ADVANCES
|
RESEARCH ARTICLE
3 of 10
despite the continuous oxygen input, indicating that
Bt
was dependent
on removal of oxygen by
Kp
. At sample points 7 (0.25 mM glucose)
and 8 (0 mM glucose),
Bt
continued to grow, despite dissolved oxygen
measurements indicating oxygen concentrations above the tolerance
for
Bt
growth (Fig. 4C). This observation differed slightly from the
model, suggesting that there may be additional biological factors beyond
metabolic coupling and stoichiometric balance of carbon and oxy-
gen that can affect multistability. Imaging revealed that in the
Kp-Bt
state, bacterial aggregates were larger at lower glucose concentrations.
Furthermore, fluorescent in situ hybridization showed these aggre-
gates contained both
Kp
and
Bt
(fig. S2). We hypothesize that coag-
gregation is one potential mechanism that could extend the region
of hysteresis by providing microenvironments more favorable for
Bt
growth by further facilitating metabolic coupling between the two
species, as observed in biofilms (
6
). Other factors, such as adhesion
to the walls of the vessel, may also contribute to extending the re-
gion of hysteresis.
Gene expression analysis of CSTR samples revealed that mul-
tistability also occurs at the transcriptome level in both the commu-
nity and in individual species. Principal components analysis (PCA)
of the community-level gene expression data showed that samples
clustered on the basis of the steady state (
Kp-
only versus
Kp-Bt
) from
which they were collected (Fig. 5A). Strong clustering at the com-
munity level is expected because
Bt
is absent from the
Kp
-only state.
However, when we evaluated the gene expression profile of
Kp
(Fig. 5B),
which is present in all steady-state conditions, we also found clustering
based on the state of the community.
To further evaluate the proposed metabolic mechanism respon-
sible for MSH (Fig. 2, B and C), we compared metabolic regulation
in
Kp
in the
Kp-Bt
state and the
Kp-
only state. We used a method
from the Neilsen Lab (
42
) to collect topological information from
the genome-scale metabolic models and combine it with gene ex-
pression data to identify reporter metabolites that maximally differ
between the two states. Among the top reporter metabolites were
pyruvate, phosphoenolpyruvate, glucose, and glucose-6-phosphate
(table S3), suggesting that the phosphotransferase system (PTS), which
is involved in sugar transport, is up-regulated in the
Kp
-only state
relative to the
Kp-Bt
state (Fig. 5, C and D). In the
Kp-Bt
state, genes
involved in the
-glucoside linked substrates were up-regulated
(Fig. 5E), suggesting that
Kp
obtains some of its carbon source from
oligosaccharides. These oligosaccharides are released into the envi-
ronment by
Bt
through the breakdown of dextran by dextranase, an
extracellular endohydrolase (
43
).
Bt
uses these oligosaccharides by
hydrolyzing them using glucan-1,3-
-glucosidases. As expected, both
dextranase (
dexA
) and glucan-1,3-
-glucosidase (
gaa
) were found
to be highly expressed in
Bt
in the
Kp-Bt
state (Fig. 5E).
A
BC
Fig. 2. A multistable model system consisting of
Kp
, a facultative anaerobe, and
Bt
, an anaerobe, that is relevant to the human gut microbiome.
(
A
) Dynamic
equations describing the model system can be solved with dFBA using each species’ genome-scale metabolic model. (
B
) In the
Kp-
only state,
Bt
does not grow and
Kp
uses external sugars and short-chain fatty acids. (
C
) In the
Kp-Bt
state,
Bt
can grow and break down complex polysaccharides into simple sugars and short-chain fatty acids,
which
Kp
can use to maintain reduced oxygen levels favorable for
Bt
growth.
on August 13, 2020
http://advances.sciencemag.org/
Downloaded from
Khazaei
et al
.,
Sci. Adv.
2020;
6
: eaba0353 12 August 2020
SCIENCE ADVANCES
|
RESEARCH ARTICLE
4 of 10
Our analysis (Fig. 5F) also suggested an up-regulation of acetate
utilization by
Kp
in the
Kp-Bt
state as inferred from the up-regulation
of acetate permease. In addition,
Kp
genes involved in lactate utili-
zation were up-regulated in the
Kp-Bt
state. Upon oxygen exposure,
Bt
is known to produce lactate (
44
). A pilot experiment showed that
the addition of a mixture of lactate and acetate can cause a direct
“jump” from the
Kp-
only state to the
Kp-Bt
state (fig. S5), emphasizing
that short-chain fatty acids are involved in the metabolic coupling
between
Kp
and
Bt
in MSH.
Multistability of gene expression
extended to the anaerobic metabolic pathway for propanediol utili-
zation in
Kp
(Fig. 5G) (
45
). We thus infer that a subpopulation of
Kp
was undergoing anaerobic metabolism in samples 4 and 5 (of the
Kp-Bt
state), where the dissolved oxygen concentrations in the reactor
were lowest (Fig. 4C). Overall, these results were consistent with
the basic mechanism for MSH (Fig. 2, B and C) and reveal that
MSH extends to the expression of genes and pathways involved in
metabolic coupling between the species.
DISCUSSION
In this work, we used genome-scale mathematical modeling, biore-
actor experiments, transcriptomics, and dynamical systems theory
to show that MSH is a mechanism that can describe shifts and per-
sistence of a two-member model microbiome aerobe-anaerobe com-
munity under seemingly paradoxical conditions (e.g., oxygen-exposed
environments). We further identified key metabolic pathways involved
AB
C
Sample point
S2S1
S3
S4
S5
S6
S7
S8
Input glucose concentration (mM)
Species concentration (millions/ml)
1500
1000
500
0
Dissolved ox
ygen
pH
Input glucose concentration
B. thetaiotaomicron
K. pneumoniae
S2S1
S3
S4
S5
S6
S7
S8
Sample point
pH
10
8
6
4
2
20
10
15
5
0
Multistable
Kp
-only state
Monostable
Kp-Bt
state
Multistable
Kp-Bt
state
Dissolved ox
ygen
concentration (uM)
0
2
4
6
Input glucose concentration (mM)
To tal cell concentration (millions/ml)
1500
1200
2
13
4
5
6
900
600
300
0
0
S4
S3
S2
S1
S5
S6
S7
S8
Kp-Bt
state
Kp
-only state
Fig. 4. MSH of
Kp
and
Bt
community in a CSTR.
(
A
) Total cell concentrations collected at the eight different steady-state sample points (S1 to S8) from the CSTR mea-
sured by qPCR. (
B
) Cell concentrations for each individual species in the community measured by qPCR. (
C
) pH and dissolved oxygen concentrations measured in the
CSTR for each sample point. Error bars are SD of at least three replicates collected at steady state (separated by >1 residence time of 5 hours) from the CSTR for each of
the eight steady-state glucose conditions. See fig. S1 for bioreactor workflow.
A
B
Input glucose concentration (mM)
Input glucose concentration (mM)
Cell conc. (millions/ml)
Cell conc. (millions/ml)
600
400
200
0
0
1
23
4
5
6
7
K. pneumoniae
300
200
0
01
23
45
67
100
B. thetaiotaomicron
B. thetaiotaomicron
Cell conc. (millions/ml)
Cell conc. (millions/ml)
800
600
400
200
0
Input oxygen flow rate (ml/min)
01234567
01234567
300
200
100
0
K. pneumoniae
Input oxygen flow rate (ml/min)
C
Input oxygen flow rate (ml/min)
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
11
.5
22
.5
33
.5
44
.5
55
.5 6
Input glucose concentration (mM)
Monostable
Kp
-
Bt
states
Monostable
Kp
-only states
Multistable states (
Kp
-only or
Kp
-
Bt
)
K. pneaumoniae
in the
Kp
-only state
B. thetaiotaomicron
in the
Kp
-only state
Kp
-
Bt
state
Fig. 3. Simulations illustrating MSH in the microbial community with respect to environmental perturbations.
Cell concentrations as a factor of (
A
) glucose
concentration variations in the input feed under constant input oxygen flow rate (1.7 ml/min) and (
B
) input oxygen flow variations under constant glucose concentrations
(3 mM) in the input feed. Each point represents the steady-state concentration for the given species in the community after a 50-hour simulation. (
C
) Regions of stability
as a function of glucose concentrations in the input feed and oxygen flow rates into the reactor. In regions of multistability (circles), the community can exist in either a
Kp
-only state or a
Kp-Bt
(aerobe-anaerobe) state under the same conditions. In regions of monostability (triangles), it is only possible for one state to exist for the given set
of parameters; down-pointed triangles represent monostable regions where
Kp
-only state exists, and up-pointed triangles represent monostable regions where only the
Kp-Bt
state can exist.
on August 13, 2020
http://advances.sciencemag.org/
Downloaded from
Khazaei
et al
.,
Sci. Adv.
2020;
6
: eaba0353 12 August 2020
SCIENCE ADVANCES
|
RESEARCH ARTICLE
5 of 10
in MSH in the
Kp-Bt
system. Future gene knockout studies would
further confirm the critical metabolic pathways responsible for MSH.
We demonstrated that altering the balance between carbon and ox-
ygen fluxes within the system, by changing input glucose levels, leads
to a community shift from the
Kp
-only state to the
Kp-Bt
state.
Follow-up experiments exploring the manipulation of input oxygen levels
would be needed to quantitatively evaluate the role of input oxygen
levels in MSH; in particular, the model-predicted ability to switch
the system back to the
Kp
-only state. A limitation of this study is the
long time frame of the CSTR experiment (fig. S1) and that the data
reported come from a single run; however, a shorter pilot experiment
(fig. S9) demonstrated similar dynamics. More broadly, identifying and
interpreting MSH in human microbiomes and microbiome-associated
diseases would require carefully designed longitudinal measurements
and models that take into account the full complexity of microbi-
omes, their spatial structure, and host responses. If MSH is found,
then it would have profound conceptual impact. To understand and
control microbial communities without MSH, one currently relies on
a well-established conceptual connection between correlation, causation,
and control. Consider points S1 to S3 (Fig. 4). The levels of
Kp
cor-
relate with the input glucose concentration; from a known input glu-
cose concentration, one can infer a steady-state
Kp
concentration
and vice versa. Input glucose concentration is the causal factor, and
therefore, it can be used to control the steady-state levels of
Kp
. If
Sample point
Gene expression (TPM)
Sample point
5000
3000
2000
S1
Gene expression (TPM)
Sample point
Sample point
Sample point
Gene expression (TPM)
Glucan 1,3-
α
-glucosidase
(
gaa
; BT3086)
Dextranase (
dexA
; BT3087)
K. pneumoniae
B. thetaiotaomicron
4000
1000
5000
4000
3000
2000
1000
0
S2 S3 S4 S5 S6 S7 S8
S1 S2 S3 S4 S5 S6 S7 S8
EF
0
G
Gene expression (TPM)
3000
2000
1000
0
K. pneumoniae
B. thetaiotaomicron
Acetate kinase (BT3693)
L
-Lactate permease (BT1453)
S1 S2 S3 S4 S5 S6 S7 S8
Gene expression (TPM)
400
300
200
100
0
S2 S3 S4 S5 S6 S7 S8
S1
Propanediol diffusion facilitator (KP03202)
Polyhedral bodies (KPN03204)
Polyhedral bodies (KPN_03203)
Dehydratase, medium subunit (KPN03206)
Dehydratase, small subunit (KPN03207)
Polyhedral bodies (KPN03211
)
Polyhedral bodies (KPN03208)
S2 S3 S4 S5 S6 S7 S8
S1
0
20
40
60
80
100
A
PC1: 45% variance
PC1: 82% variance
PC2: 25% variance
PC2: 4% variance
S2 (1 mM)
S3 (2 mM)
S4 (5 mM)
S5 (2 mM)
S6 (1 mM)
S7 (0.25 mM)
S8 (0 mM)
0
-
10
-
20
-
10
10
Kp
-only state
Community transcriptome
K. pneumoniae
transcriptome
20
0
-
20
S1 (0.25 mM)
0
10
20
-
100
100
0
-
40
-
60
Kp
-only state
B
Glucose
Glucose-6-
P
Mannose-6-
P
Pyruvate
Phosphoenol-
pyruvate
C
ptsG
KPN01099
crr
KPN02764
ptsI
KPN02763
Periplasm
glyC/glyB
KPN04086
manZ
KPN02335
manX
KPN02333
P
P
Mannose
ptsH
KPN02762
P
Gene expression (TPM)
Gene expression (TPM)
Gene expression (TPM)
ptsI
ptsH
ptsG
crr
manZ
manX
12
3451
2
345
12
345
0
00
0
0
0
1000
2000
3000
4000
2000
4000
6000
3000
6000
9000
12000
P
D
Input glucose
concentration (mM)
Input glucose
concentration (mM)
Input glucose
concentration (mM)
PTS family enzyme IIC
(KPN04086)
Glycoside hydrolase, family 4
(KPN04085)
L
-Lactate permease (KPN03947)
L
-Lactate dehydrogenase
(KPN03949)
Acetate permease
(KPN04476)
Kp-Bt
state
Kp-Bt
state
Fig. 5. Gene expression analysis of CSTR steady-state samples.
(
A
) PCA of the community transcriptome; each point represents the combined transcriptome of
Kp
and
Bt
for each sample (S1 to S8). (
B
) PCA of the
Kp
transcriptome. (
C
) The most differentially regulated pathway between the
Kp-
only and the
Kp-Bt
states is the PTS.
The
gray
box indicates the up-regulated gene; white boxes are down-regulated genes. (
D
) PTS genes down-regulated in the
Kp-Bt
state in
Kp
. Solid lines represent the
Kp-
only
state, and dashed lines represent the
Kp-Bt
state. (
E
) Gene expression, in transcripts per million (TPM), of oligosaccharide uptake in
Kp
and dextran metabolism to oligo-
saccharides in
Bt
for each steady-state sample point. (
F
) Expression of genes involved in acetate and lactate utilization in
Kp
and acetate and lactate production in
Bt
for
each CSTR sample. (
G
) Expression of the propanediol utilization pathway in
Kp.
(E to G) Unshaded regions are the
Kp-
only state; the gray shaded region is the
Kp-Bt
state.
Error bars are SD of three replicates collected (separated by >1 residence time) from the CSTR for each of the eight steady-state glucose conditions.
on August 13, 2020
http://advances.sciencemag.org/
Downloaded from
Khazaei
et al
.,
Sci. Adv.
2020;
6
: eaba0353 12 August 2020
SCIENCE ADVANCES
|
RESEARCH ARTICLE
6 of 10
MSH is identified in microbiomes, then it would break this familiar
conceptual connection between causation and correlation. Consider
the region of hysteresis (points S1 to S3 and S5 to S7; Fig. 4). The
observed steady-state levels of
Kp
no longer correlate with the input
glucose concentration. At 2 mM input glucose, the system could be
in either the
Kp
-only state S3 or the
Kp-Bt
state S5. At ~650 × 10
6
colony-forming units/ml of
Kp
, the input glucose levels could be
either 0.25 or 2 mM.
Although there is no correlation between species’
abundance and input glucose concentration (in other words, knowing
glucose concentration is not sufficient for predicting species abun-
dance; instead, a system’s history must also be known), input glucose
concentration remains the causal factor. Furthermore, under MSH,
establishing causation is insufficient for achieving control: Although
input glucose concentration is the causal factor responsible for changes
in the community state, it cannot be used to fully control the com-
munity (i.e., one cannot use changes in glucose inputs to revert the
Kp-Bt
state back to the
Kp
-only state). Alternative control strategies
(e.g., changes in oxygen levels or disruption of metabolic coupling),
derived from appropriate models, would need to be deployed under
MSH.
Therefore, recognizing whether and when MSH exists in human
microbiomes will be critical for interpreting correlation and causation
and for designing therapeutic control strategies that can steer mi-
crobial communities to desirable states.
MATERIALS AND
METHODS
Model development
For the computational simulations, we used the DMMM framework
(
38
), which is an extension of dFBA applied to microbial communities.
Briefly, the DMMM framework has two components. In component 1,
external differential equations describe mass balances for species
and metabolite concentrations in the CSTR (shown in Fig. 2A and
described in the Supplementary Materials). Unlike the traditional
method for solving differential equations in a bacterial system, we
do not assume that parameters such as growth rates and metabolic
flux rates are constant. Instead, we allow the parameters to be dynamic
because we are studying a system with potentially rich dynamics. To
find the values for these dynamic parameters, we use FBA (compo-
nent 2) to solve for the parameter values at every time step of the
simulated time period. Component 2 includes the genome-scale
models for each species. These models are used to perform FBA at
every time point to obtain updated parameters for the differential
equations in component 1.
Component 1
The system is described as a CSTR with the following mathematical
formulation
d
X
i
dt
=
i
X
i
F
out
X
i
V
(1)
dS
j
dt
=
i
v
i
j
X
i
+
F
in
S
feed
j
F
out
S
j
V
(2)
dS
oxygen
dt
=
i
v
i
oxygen
X
i
+
K
L
a
(
S
*
S
oxygen
)
(3)
Here,
V
is the volume of the reactor (constant), and
X
i
is the
biomass (g/liter) of the
i
th microbial species.
S
j
is the concentration
(mM) of the
j
th metabolite,
F
in
is the rate of flow (liter/hour) into the
reactor,
F
out
is the rate of flow (liter/hour) out of the reactor,
S
feed
j
is
the concentration of the
j
th metabolite in the feed stream,
i
(hour
−1
)
is the growth rate of the
i
th microbial species,
v
i
j
is the metabolic
flux of the
j
th substrate in the
i
th microbial species,
K
L
a
is the volu-
metric oxygen transfer coefficient, and
S
*
is the dissolved oxygen
saturation concentration.
In continuous culture, substrate utilization can deviate from di-
auxic growth (as typically observed in batch culture) and cosubstrate
utilization is possible (
46
). The set of differential equations are solved
using the following analytical approximation
X
f
=
X
i
,0
e
(
i
F
in
_
V
0
)
T
(4)
S
f
j
=
S
0
j
+
i
[
v
i
j
V
0
μ
i
V
0
F
in
(
X
i
,0
e
(
μ
i
F
in
_
V
0
)
Δ
T
X
i
,0
)
]
+
F
in
(
S
feed
j
S
j
)
V
0
Δ
T
(5)
S
f
oxygen
=
S
0
oxygen
+
i
[
v
i
oxygen
V
0
μ
i
V
0
F
in
(
X
i
,0
e
(
μ
i
F
in
_
V
0
)
Δ
T
X
i
,0
)
]
+
K
L
a
(
S
*
S
oxygen
)
Δ
T
(6)
At the beginning of every time step (∆
T
), the parameters
i
and
v
i
j
are calculated using FBA from genome-scale models (component
2) and fed back into Eqs. 4, 5, and 6. This process is repeated for all
time intervals in the simulated time period.
Component 2
Genome-scale metabolic models are used to establish genotype-
phenotype relationships and capture the metabolic capabilities of
each model organism. Furthermore, these models allow us to inte-
grate metabolic network topology information with RNA sequenc-
ing data for transcriptomic analysis of the CSTR experiments
(described further in the “RNA sequencing and analysis” section).
We used the published iYL1228 model of
Kp
MGH-78578 (
37
) and
the published iAH991 model of
Bt
VPI-5482 (
36
). Because these
models are well validated (and curated), we added the minimum
number of parameters that allow for integration of these models to
the community dFBA framework. Our changes include adding a
pathway for dextran uptake and hydrolysis to glucose in the
Bt
iAH991 model. The pathway lumps hydrolysis of dextran to glucose
into a single reaction. In this lumped reaction, we assume that 50%
of the glucose produced from dextran by
Bt
can be released into the
environment for shared use. For the purpose of the simulations, dex-
tran is assumed to be 100 glucose units. The genome-scale models are
solved separately for each species by FBA (
47
) at each time point
max
c
T
̄
v
i
(7)
s
.
t
.
A
i
̄
v
i
=
0
̄
v
i
,lb
<
̄
v
i
<
̄
v
i
,ub
where
c
is the cost vector,
̄
v
is the vector of fluxes, and
A
is the ma-
trix of mass balance stoichiometries. The optimization criterion is
biomass growth rate (for each species). For the bounds for the fluxes,
we used the values in the curated, published genome-scale models
on August 13, 2020
http://advances.sciencemag.org/
Downloaded from
Khazaei
et al
.,
Sci. Adv.
2020;
6
: eaba0353 12 August 2020
SCIENCE ADVANCES
|
RESEARCH ARTICLE
7 of 10
(
Kp
iYL1228 and
Bt
iAH99). The uptake fluxes explicitly modeled
in component 1 (dextran, glucose, acetate, and oxygen) are bounded
by Michaelis-Menten kinetics
v
i
,ub
j
=
v
i
j
,max
S
j
K
m
+
S
j
(8)
The values for
v
i
j
,max
and
K
m
for some of the metabolites in the
model were estimated from batch experiments. Batch culture ex-
periments were carried out in a 96-well flat-bottom plate. Overnight
cultures grown anaerobically in minimal medium supplemented
with either 0.5% (w/v) dextran or 0.5% (w/v) glucose were diluted
1:20 (for
Bt
) and 1:100 (for
Kp
) and outgrown to mid-log phase. The
cultures were then pelleted and resuspended at optical density (OD)
1 (for
Bt
) and OD 0.1 (for
Kp
) in carbon-source-free minimal medium.
We added 10
l of cells to 200
l of minimal medium containing
various concentrations [0.125 to 0.5% (w/v)] of the carbon source.
The plate was incubated at 37°C, and OD
600
was measured every
10 min. For batch cultures, Monod growth kinetics was assumed
dX
dt
=
X
max
S
S
+
K
(9)
dS
dt
=
v
max
S
S
+
K
(10)
where
max
is the maximum growth rate the given bacteria can
achieve on a carbon source when it is not resource limited. Growth
data from replicate wells of multiple concentrations of carbon source
were fitted simultaneously using Bayesian parameter estimation
implemented with Markov chain Monte Carlo (
48
). Individual
growth curves were allowed to have distinct initial cell concentrations
and background values, with other parameters held constant. The
fitted parameters are presented in table S1 and fig. S3. For all other
metabolites captured in the differential equations, the
K
m
and
v
max
values are assumed to be the same (
v
max
of 10 mmol/gCDW·hour
and
K
m
of 0.01 mM) on the basis of literature for
Escherichia coli
(
49
,
50
).
Values for parameters and initial conditions used in the model
are presented in table S2. Initial conditions are chosen to represent
the experimental setup, whereby we first establish a steady state for
Kp
in the CSTR before inoculating
Bt
. Therefore, in the models, we
start with a higher concentration of
Kp
than
Bt
. The initial conditions
for
Kp
in the reactor is arbitrarily chosen to be the experimentally
measured monoculture steady-state concentration of
Kp
at an input
glucose concentration of 0.25 mM.
The initial conditions for
Bt
in
the reactor are 0.0015 g/liter, which is equivalent to addition of 1 ml of
OD 1
Bt
into the reactor, as done experimentally. For most steady-
state conditions, glucose is limiting, and therefore, the initial condi-
tions for glucose concentration in the reactor are chosen to be 0 mM.
We note that although an individual FBA simulation (component
1) is linear, a dynamic FBA (the combination of component 1 and
component 2) is not linear. The behavior of the cells will shift mark-
edly (a nonlinear response) when substrate levels in the environment
cross certain thresholds. The two biggest sources of nonlinearity in
our model system are binary growth/no-growth behavior of
Bt
in
the presence or absence of O
2
and the binary capacity of
Kp
to use
carbon from dextran (through its breakdown into oligosaccharides
by
Bt
) in the presence or absence of
Bt
. To model the binary growth/
no-growth behavior of
Bt
in the presence of O
2
, we included a con-
ditional operator before the FBA simulations for each time step: If the
O
2
levels are above 350 nM, then the growth rate of
Bt
is set to zero
and the FBA is not run for
Bt
, whereas if O
2
levels are lower than
350 nM O
2
, then the FBA simulation for
Bt
is allowed to proceed. The
value for the O
2
growth threshold for the model is arbitrarily chosen
from the concentration range reported for the aerotolerant
Bt
(
51
).
To computationally identify the regions of stability with respect
to glucose and oxygen (Fig. 3C and fig. S4), we varied oxygen input
flow rates at constant input glucose concentration for each glucose
condition examined. We evaluated 11 glucose conditions ranging from
1 to 6 mM.
For each given glucose input concentration, we started
with oxygen at an input flow rate of 6 ml/min and ran the simula
-
tion for 50 hours to ensure that the system reached a steady state.
We then decreased the oxygen input by intervals of 0.5 ml/min down
to 0.5 ml/min, for each oxygen condition, ensuring that the system
reached a steady state [we refer to the oxygen variations at 6 to
0.5 ml/min for a given constant glucose input concentration as the
“forward simulations”]. The concentration of oxygen input at which
Bt
starts to grow is identified as the tipping point to the monostable
Kp-Bt
state. After running the oxygen simulation at 0.5ml/min,
we increased the concentration back to 6 ml/min at intervals of
0.5 ml/min [we refer to the oxygen variations at 0.5 to 6 ml/min as
the “reverse simulations”). The concentration of oxygen at which
Bt
can no longer grow and gets washed out is identified as tipping
point to the monostable
Kp-
only state. The region between these
two tipping points to the monostable
Kp-Bt
state in the forward
simulations and the monostable
Kp
-only state in the reverse simu-
lations is identified as the region of bistability. The colors in fig. S4
represent the steady-state concentration of
Kp
in the reverse simu-
lations divided by the steady-state concentration of
Kp
in the
forward simulations. In regions of monostability, the concentration
of
Kp
is similar in both the forward and reverse simulations and
therefore has a value of approximately 1.
Continuous culture of
Kp
and
Bt
Continuous culture experiments were carried out in a 500-ml
bioreactor (miniBio Applikon Biotechnology, Delft, Netherlands)
with a total culture volume of 200 ml. Minimal medium [KH
2
PO
4
(3.85 g/liter), K
2
HPO
4
(12.48 g/liter), (NH
4
)
2
SO
4
(1.125 g/liter),
1× methyl methanesulfonate (20× methyl methanesulfonate: NaCl,
17.6 g/liter; CaCl
2
, 0.4 g/liter; MgCl
2
× 6H
2
O, 0.4 g/liter; MnCl
2
× 4H
2
O,
0.2 g/liter; and CoCl
2
× 6H
2
O, 0.2 g/liter), Wolfe’s mineral solution
(10 ml/liter) (
52
), Wolfe’s vitamin solution (10 ml/liter) (
52
), 4.17
M
FeSO
4
× 7H
2
0, 0.25 mM cysteine, 1
M menadione, 2
M resazurin,
dextran (1 g/liter; Sigma, D5376; average molecular weight, 1.5 × 10
6
to 2 × 10
6
), and glucose at varying concentrations] was purged with
100% N
2,
stored under anaerobic conditions before use, and main
-
tained under N
2
during operation of the CSTR.
We calibrated the
dissolved oxygen probe by aerating the reactor with CO
2
(5 ml/min)
and N
2
(45 ml/min). The stable measurement without O
2
input was
taken to be 0
M dissolved oxygen. We then sparged the reactor with
O
2
(1.7 ml/min), CO
2
(5 ml/min), and N
2
(43.3 ml/min). This stable
measurement was taken to be 30.714
M O
2
[calculated assuming
dissolved oxygen (6.056 mg/liter) at 1 atm and 37°C].
For the experiments, the bioreactor was sparged with total gas at
50 ml/min [O
2
(1.7 ml/min), CO
2
(5 ml/min), and balance of N
2
]
and agitated with two six-bladed Rushton turbines operated at
750 rpm. Temperature was maintained at 37°C, and a residence time
of 5 hours (input and output flow rates of 40 ml/hour) was used for
on August 13, 2020
http://advances.sciencemag.org/
Downloaded from