Page
1
of
20
1
2
3
Metabolic multi-stability and hy
steresis in a model aerobe-anae
robe
4
microbiome community
5
6
Tahmineh Khazaei,
1
Rory L. Williams,
1
Said R. Bogatyrev,
1
John C. Doyle,
1,2
7
Christopher S. Henry,
3
Rustem F. Ismagilov
1,4*
8
9
1
Division of Biology and Biologica
l Engineering, California Inst
itute of Technology
10
2
Division of Engineering and Applied S
ciences, California Instit
ute of Technology
11
3
Data Science and Learning Divisio
n, Argonne National Laboratory
12
4
Division of Chemistry and Chemical Engineering, California Inst
itute of Technology
13
*Correspondence to:
rustem.admin@caltech.edu
14
15
16
17
18
19
Abstract
20
21
Changes in the composition of the human microbiome are associat
ed with health and disease. Some
22
microbiome states persist in seemingly unfavorable conditions,
e.g., the proliferation of aerobe-
23
anaerobe communities in oxygen-exposed environments in wounds o
r small intestinal bacterial
24
overgrowth. However, it remains
unclear how different stable mi
crobiome states can exist under
25
the same conditions, or why some states persist under seemingly
unfavorable conditions. Here,
26
using two microbes relevant to the human microbiome, we combine
genome-scale mathematical
27
modeling, bioreactor experiments, transcriptomics, and dynamica
l systems theory, to show that
28
multi-stability and hysteresis (MSH) is a mechanism that can de
scribe the shift from an aerobe-
29
dominated state to a resilient, paradoxically persistent aerobe
-anaerobe state. We examine the
30
impact of changing oxygen and nutrient regimes and identify fac
tors, including changes in
31
metabolism and gene expression, that lead to MSH. When analyzin
g the transitions between the
32
two states in this system, the familiar conceptual connection b
etween causation and correlation is
33
broken and MSH must be used to interpret the dynamics. Using MS
H to analyze microbiome
34
dynamics will improve our conceptual understanding of the stabi
lity of microbiome states and the
35
transitions among microbiome states.
36
37
38
39
One sentence summary
40
Multi-stability and hysteresis (M
SH) is a potential mechanism t
o describe shifts to and persistence
41
of aerobe-anaerobe communities in the microbiome.
42
43
44
45
46
47
48
49
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
2
of
20
Introduction
50
51
Recent evidence shows that changes in the species composition a
nd abundance of the human
52
microbiome can be associated with health and disease (
1-3
). Understanding the mechanisms that
53
cause compositional shifts in healthy microbiomes, which otherw
ise can be remarkably stable, is
54
challenging due to the inherent c
omplexity of these ecosystems. A perplexing feature of some of
55
these disturbed ecosystems is the persistence of a new microbio
me state, even in seemingly
56
unfavorable conditions. For example, in small intestinal bacter
ial overgrowth (SIBO), strict
57
anaerobes that are typically f
ound only in the colon become pro
minent in the small intestine and,
58
paradoxically, persist in this
environment exposed to oxygen fl
ux from the tissue (
4, 5
). Similarly,
59
in periodontal diseases (
6
) and in wound infections, anaer
obes proliferate in oxygen-expo
sed
60
environments.
61
62
One potential mechanism to explain microbiome shifts and their
persistence is multi-
63
stability (
7-13
), the concept that several steady states can exist for an identical set of system
64
parameters (Fig. 1).
65
66
Fig. 1. Cartoon explanation of multi-stability without hysteres
is (A) and with hysteresis (B) using a seesaw and
67
child analogy. (A)
Consider an idealized seesaw. When a child stands on the left s
ide of the seesaw (positions 1, 2), the
68
seesaw is tilted to the left (the “left” state). When the chil
d walks past the mi
dpoint (positions 3, 4), this idealized
69
seesaw would tilt to the right (the “right” state). The seesaw
is a multi-stable system that can be in either of the two
70
states depending on the position of the child. Moreover, the st
ate of this seesaw is correlated with the position of the
71
child.
(B)
Hysteresis arises if the seesaw
is rusted and “sticky.” As the
child walks through positions 1 and 2 past the
72
midpoint to position 3, the seesa
w remains stuck in the left-st
ate; the child has to walk much further, to position 4,
73
to switch it to the right state. As the child walks back to th
e left, through positions 5, 6, and past the midpoint to positi
on
74
7, the seesaw remains stuck in th
e right state and it only swit
ches back to the left state when the child reaches position
75
8. In the region of positions 2, 3, 6, and 7, the system exhib
its multi-stability and hysteresis. Under MSH there is no
76
correlation between the causal f
actor (position of the child) a
nd the observed state of the seesaw. Furthermore, the
77
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
3
of
20
causal factor (position of the child) cannot be used to control
the state of the seesaw within the region of MSH. The
78
system must escape the region of
MSH to switch the state of the
seesaw. Applying the cartoo
n to this manuscript, the
79
states of the model microbiome
become analogous to the state
s of the seesaw, and the changes in the various inputs,
80
such as carbon and oxygen, become analogous to the changes in t
he child’s position.
81
82
Multi-stable systems have been described in the context of ecos
ystems (
14-17
) and gene-regulatory
83
networks (
18-20
). Now, with the expanding characterization of the microbiome, there are signs that
84
multi-stability may also exist in these communities (
21-27
). For example, compositional changes
85
in gut microbiota are implicated in inflammatory bowel disease
(
28
) and obesity (
29
). Bimodal
86
species abundance (i.e., when a microbial species is present at
either high or low levels) has been
87
interpreted as multi-stability (
30
); however, as discussed by Gonze et al., bimodality is insuffi
cient
88
to prove multi-stability (
7, 31
). Some multi-stable systems can additionally exhibit hysteresi
s,
89
where in response to a perturbation, a system gets “stuck” in a
new steady state and the former state
90
cannot be regained by simply reversing the perturbation (
7
). The presence of hysteresis could be
91
hypothesized from studies of the microbiome (
32
). For example, antibiotic exposures can change
92
the microbiome composition, and ha
ve lasting effects even after
removal of the antibiotic (
33, 34
).
93
However, it has not been rigorously tested whether multi-stabil
ity and hysteresis (MSH) can arise
94
in a microbiome-relevant community and by what mechanism.
95
Here, we investigate MSH in a minimally “complex” two-species s
ystem to represent the
96
paradoxical aerobe–anaerobe microbiome communities that persist in oxygen-exposed
97
environments. We used two organisms prevalent in SIBO (
35
): the anaerobe
Bacteroides
98
thetaiotaomicron
(
Bt
) that breaks down complex carbohydrates (e.g. dextran) into si
mple sugars
99
and short chain fatty acids (
36
), and the facultative anaerobe (hereafter referred to as an
aerobe)
100
Klebsiella pneumoniae
(
Kp
) capable of consuming oxygen, simple sugars and short chain fa
tty
101
acids, and performing anaerobic re
spiration in the absence of o
xygen (
37
).
102
103
Results
104
105
To simulate how the interplay between environmental perturbatio
ns and inter-species metabolic
106
interactions could lead to multi-stability, we first built a ma
thematical model. We used the dynamic
107
multi-species metabolic modeling (DMMM) framework (
38
) to model a community of
Kp
and
Bt
108
in a continuously stirred tank rea
ctor (CSTR) (Fig. 2A) with co
ntinuous input flows of dextran
109
minimal media and varying input glucose or varying input oxygen
levels (depending on the
110
simulation). The outflow rate is
equal to the inflow rate to ma
intain a constant reactor volume. The
111
DMMM framework uses dynamic flux balance analysis (dFBA) (
34
), which allows us to capture
112
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
4
of
20
temporal changes in intracellular flux rates (using the genome-
scale metabolic model for each
113
species), extracellular metabol
ite concentrations, and species
concentrations.
114
115
116
Fig. 2. A multi-stable model system consisting of
Klebsiella pneumoniae (Kp)
, a facultative anaerobe, and
117
Bacteroides thetaiotaomicron (Bt)
, an anaerobe, that is relevant
to the human gut microbiome.
(
A
) Dynamic
118
equations describing the model system can be solved with dynamic flux-balance analysis
utilizing each species’
119
genome-scale metabolic model. (
B
) In the
Kp-
only state,
Bt
does not grow and
Kp
utilizes external sugars and short
120
chain fatty acids. (
C
) In the
Kp–Bt
state,
Bt
can grow and break down complex
polysaccharides into simple su
gars and
121
short chain fatty acids, which
Kp
can utilize to maintain reduced oxygen levels favorable for
Bt
growth.
122
123
Next, to computationally test whe
ther a nutrient
perturbation c
ould lead to a cha
nge in community
124
state, we altered glucose input concentrations (Fig. 3A), while
keeping constant all other system
125
parameters, including continuous
oxygen input and continuous de
xtran input. The model predicted
126
that for glucose concentrations
of 0.25–3 mM in the input feed
(at a constant flow rate of 0.7
127
mL/min for all conditions), the output state consisted solely o
f
Kp
, which we refer to as the
Kp
-
128
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
5
of
20
only state (Fig. 2B)
.
Stoichiometrically, at these glucose concentrations oxygen was
not completely
129
consumed, thus the environment was unfavorable for
Bt
growth. However, when we increased
130
glucose input concentration to
3.25 mM, we observed a shift to
a new steady state (Fig. 3A). At this
131
“tipping point,” the environment became sufficiently anaerobic
to support the growth of
Bt
. We
132
refer to this second distinct steady state as the
Kp–Bt
(aerobe–anaerobe) state (Fig. 2C). In the
Kp–
133
Bt
state, the
Kp
population is no longer carbon limited due to the additional c
arbon sources
134
generated from the metabolism of dextran by
Bt
. Thus,
Kp
can now consume all of the available
135
oxygen to oxidize both glucose and the additional carbon source
s, resulting in anaerobic conditions.
136
Surprisingly, this
Kp–Bt
state persisted even when we sy
stematically rev
ersed the input of glucose
137
below 3.25 mM, even to 0 mM. Thus, this system shows hysteresis
and multi-stability: under
138
identical input conditions of glucose and oxygen, the system ca
n be in either of the two possible
139
states. We then identified tipping points for population shifts
in response to input oxygen
140
variations—with glucose kept constant (Fig. 3B
). In addition to glucose input, all other parameters,
141
including continuous dextran input, are held constant.
We considered oxygen as a parameter because
142
in host settings oxygen availability can be affected by perfusi
on rate, blood flow rate, immune
143
consumption, etc. We found that we could return the system to t
he
Kp
-only state by increasing
144
oxygen levels, a state switch that
was not possible by manipula
ting glucose concentration alone.
145
Finally, we simulated changes in both glucose and oxygen levels and characterized the landscape
146
of multi-stability and mono-stability in the model microbial co
mmunity (Fig. 3C). These simulation
147
results illustrate that even a minimal model of microbiome with
co-
dependence (
39, 40
) can
148
demonstrate dramatic MSH.
149
150
Fig. 3. Simulations illustrating
multi-stability and hysteresis
in the microbial community with respect to
151
environmental p
erturbations.
Cell concentrations as a factor of (
A
) glucose-concentration variations in the input feed
152
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
6
of
20
under constant input oxygen flow rate (1.7 mL/min), and (
B
) input oxygen flow variations under constant glucose
153
concentrations (3 mM) in the input feed. Each point represents
the steady-state concentration for the given species in
154
the community after a 50-h simulation. (
C
) Regions of stability as a function of glucose concentrations
in the input
155
feed and oxygen flow rates into the reactor. In regions of mult
i-stability (circles) the community can exist in either a
156
Kp
-only state or a
Kp–Bt
(aerobe–anaerobe) state under th
e same conditions. In regions of mono-stability (triangles) it
157
is only possible for one state to exist for the given set of pa
rameters; down-pointed triangles represent mono-stable
158
regions where
Kp
-only state exists and up-pointed triangles represent mono-stab
le regions where only the
Kp–Bt
state
159
can exist.
160
161
We next tested these computationa
l predictions experimentally i
n a CSTR (200 mL culture
162
volume), and further explored the
metabolic factors behind the
dynamics of this aerobe–anaerobe
163
community. We varied glucose co
ncentrations (while keeping all
other parameters constant,
164
including continuous dextran and oxygen)
and measured the steady state output composition of the
165
microbial community by qPCR (and digital PCR; Fig. S6). Oxygen
was sparged into the reactor at
166
3.4% of the gas feed (50 mL/min total gas feed) and kept consta
nt for all conditions. For each
167
steady-state condition, we collected three CSTR samples separat
ed by at least one residence time
168
(5 h) (Fig. S1 contains the experimental workflow).
169
As predicted by the mathematical models, we observed both multi
-stability and hysteresis
170
(Fig. 4A) experimentally. At 0.25 mM, 1 mM, and 2 mM glucose co
ncentrations, the steady-state
171
community consisted only of
Kp
;
Bt
was washed out under these conditions (Fig. 4B). The
172
dissolved-oxygen measurements (Fig. 4C) confirmed that oxygen w
as not limiting under the
173
selected parameter conditions, re
sulting in an aerobic environm
ent unsuitable for
Bt
growth. As in
174
the simulations, at 5 mM glucose, a new distinct steady state w
as reached where
Bt
grew in the
175
presence of
Kp.
Although there was continuous oxygen flux into the reactor, th
e concentration of
176
dissolved oxygen measured in the reactor was near zero. Next, t
o test for hysteresis, we reduced the
177
glucose input back down to 2 mM, 1 mM, 0.25 mM, and 0 mM and fo
und that the aerobe–anaerobe
178
state persisted. The persistence of the
Kp–Bt
state (instead of a return to the
Kp
-only state),
179
qualitatively confirmed model pr
edictions of hysteresis, and ve
rified that this microbial community
180
is a multi-stable system.
181
The CSTR results demonstrate metabolic coupling and co-depende
nce between these two
182
bacterial species with respect
to carbon and oxygen. At sample point 8, there is no glucose input to
183
the reactor, yet
Kp
continued to grow, indicating that
Kp
was completely dependent on
Bt
for its carbon
184
supply. At sample point 4,
Bt
started to grow, despite the continuous oxygen input, indicating that
Bt
185
was dependent on removal of oxygen by
Kp
. At sample points 7 (0.25 mM glucose) and 8 (0 mM
186
glucose),
Bt
continued to grow, despite dissolved-oxygen measurements indic
ating oxygen
187
concentrations above the tolerance for
Bt
growth (Fig. 4C). This observation differed slightly from the
188
model, suggesting that there may be additional biological facto
rs beyond metabolic coupling and
189
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
7
of
20
stoichiometric balance of car
bon and oxygen that can affect mul
ti-stability. Imaging revealed that in
190
the
Kp–Bt
state bacterial aggregates were larger at lower glucose concent
rations. Furthermore,
191
fluorescent
in situ
hybridization (FISH) showed these aggregates contained both
Kp
and
Bt
(Fig. S2).
192
We hypothesize that co-aggregation is one potential mechanism t
hat could extend the region of
193
hysteresis by providing microenvironments more favorable for
Bt
growth by further facilitating
194
metabolic coupling between the two species, as observed in biof
ilms (
6
). Other factors, such as
195
adhesion to the walls of the vesse
l may also contribute to exte
nding the region of hysteresis.
196
197
198
Fig. 4.
Multi-stability and hysteresis of
K. pneumoniae
(
Kp
) and
B. thetaiotaomicron
(
Bt
) community in a CSTR
.
199
(
A
) Total cell concentrations collected at the eight different st
eady state sample points (S1–S8) from the CSTR
200
measured by qPCR. (
B
) Cell concentrations for each individual species in the commun
ity measured by qPCR. (
C
) pH
201
and dissolved-oxygen concentrations measured in the CSTR for ea
ch sample point. Error bars are S.D. of at least three
202
replicates collected at steady st
ate (separated by >1 residence
time of 5 h) from the CSTR
for each of the eight steady-
203
state glucose conditions. See Fi
g. S1 for bioreact
or workflow.
204
205
Gene-expression analysis of CSTR samples revealed that multi-stability also occurs at the
206
transcriptome level in both the community and in individual spe
cies. Principal component analysis
207
(PCA) of the community-level gene expression data showed that s
amples clustered based on the
208
steady state (
Kp-
only vs.
Kp–Bt)
from which they were collected (Fig. 5A). Strong clustering at
the
209
community level is expected because
Bt
is absent from the
Kp
-only state. However, when we
210
evaluated the gene-expression profile of
Kp
(Fig. 5B), which is present in all steady state conditions,
211
we also found clustering based on t
he state of the community.
212
To further evaluate the proposed metabolic mechanism responsib
le for MSH (Fig. 2B,C),
213
we compared metabolic regulation in
Kp
in the
Kp–Bt
state and the
Kp-
only state.
We used a method
214
from the Neilsen lab (
41
) to collect topological informa
tion from the genome-scale meta
bolic
215
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
8
of
20
models and combine it with gene-e
xpression data to identify rep
orter metabolites that maximally
216
differ between the two states. Among the top reporter metabolit
es were pyruvate,
217
phosphoenolpyruvate, glucose, and glucose-6-phosphate (Table S3
), suggesting that the
218
phosphotransferase system (PTS), which is involved in sugar tra
nsport, is upregulated in the
Kp
–
219
only state relative to the
Kp-Bt
state (Fig. 5C–D). In the
Kp–Bt
state, genes involved in the alpha-
220
glucoside linked substrates were upregulated (Fig. 5E), suggest
ing that
Kp
obtains some of its
221
carbon source from oligosaccharides
. These oligosaccharides are
released into the environment by
222
Bt
through the breakdown of dextran by dextranase, an extracellul
ar endohydrolase (
42
).
Bt
utilizes
223
these oligosaccharides by hydrolyzing them using glucan-1,3-alp
ha-glucosidases. As expected,
224
both dextranase (
dexA
) and glucan-1,3-alpha-glucosidase (
gaa
) were found to be h
ighly expressed
225
in
Bt
in the
Kp–Bt
state (Fig. 5E).
226
Our analysis (Fig. 5F) also suggested an upregulation of aceta
te utilization by
Kp
in the
Kp–
227
Bt
state as inferred from the upregulation of acetate permease. A
dditionally,
Kp
genes involved in
228
lactate utilization were upregulated in the
Kp–Bt
state. Upon oxygen exposure,
Bt
is known to
229
produce lactate (
43
). A pilot experiment showed that the addition of a mixture of
lactate and acetate
230
can cause a direct “jump” from the
Kp-
only state to the
Kp–Bt
state (Fig. S5), emphasizing that
231
short-chain fatty acids are involved in the metabolic coupling
between
Kp
and
Bt
in MSH. Multi-
232
stability of gene expression exte
nded to the anaerobic metabolic pathway for propanediol utilization
233
in
Kp
(Fig. 5G) (
44
). We thus infer tha
t a subpopulation of
Kp
was undergoing anaerobic
234
metabolism in samples 4 and 5 (of the
Kp–Bt
state), where the dissolved
oxygen concentrations in
235
the reactor were lowest (Fig. 4C
). Overall, these results were
consistent with the basic mechanism
236
for MSH (Figure 2B–C) and reveal that MSH extends to the expres
sion of genes and pathways
237
involved in metabolic coup
ling between the species.
238
239
240
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
9
of
20
241
Fig. 5.
Gene-expression analysis of CSTR steady-state samples.
(
A
) PCA of the community transcriptome; each
242
point represents the combined transcriptome of
K. pneumoniae
(
Kp
) and
B. thetaiotaomicron
(
Bt
) for each sample (S1-
243
S8). (
B
) PCA of the
Kp
transcriptome. (
C
) The most differentially regulated pathway between the
Kp-
only and the
Kp–
244
Bt
states is the phosphotransferase
system (PTS); the grey box ind
icates the upregulated gene; white boxes are
245
downregulated genes. (
D
) PTS genes downregulated in the
Kp–Bt
state in
Kp
. Solid lines represent the
Kp-
only state
246
and dashed lines represent the
Kp–Bt
state. (
E
) Gene expression, in transcripts per million (TPM), of oligosa
ccharide
247
uptake in
Kp
and dextran metabolism to oligosaccharides in
Bt
for each steady-state sample point. (
F
) Expression of
248
genes involved in acetate an
d lactate utilization in
Kp,
and acetate and l
actate production in
Bt
for each CSTR sample.
249
(
G
) Expression of the propanediol-utilization pathway in
Kp.
(
E–G
) Unshaded regions are the
Kp-
only state; the gray
250
shaded region is the
Kp–Bt
state. Error bars are S.D. of th
ree replicates collected (separ
ated by >1 residence time) from
251
the CSTR for each of
the eight steady-state glucose conditions.
252
253
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
10
of
20
Discussion
254
255
In this work, we used genome-scale mathematical modeling, biore
actor experiments,
256
transcriptomics, and dynamical systems theory to show that MSH
is a mechanism that can describe
257
shifts and persistence of a two-member model microbiome aerobe–
anaerobe community under
258
seemingly paradoxical conditions (e.g., oxygen-exposed environm
ents). We further identified key
259
metabolic pathways involved in MSH in the
Kp-Bt
system. Future gene knock-out studies would
260
further confirm the critical met
abolic pathways responsible for
MSH.
A limitation of this study is
261
the long timeframe of the CSTR experiment (Fig. S1) and that th
e data reported come from a single
262
run; however, a shorter pilot experiment (Fig. S9) demonstrated
similar dynamics.
More broadly,
263
identifying and interpreting MSH in human microbiomes and micro
biome-associated diseases
264
would require carefully designe
d longitudinal measurements and models that take into account the
265
full complexity of microbiomes, their spatial structure, and ho
st responses. If MSH is found, it
266
would have profound conceptual impact. To understand and contr
ol microbial communities
267
without MSH, one currently relies
on a well-established concept
ual connection between correlation,
268
causation, and control. Consider points S1-S3 (Fig. 4). The le
vels of
Kp
correlate with the input
269
glucose concentration—from a known input glucose concentration,
one can infer a steady-state
Kp
270
concentration and vice versa. Inpu
t glucose concentration is th
e causal factor and therefore it can
271
be used to control the steady-state levels of
Kp
. If MSH is identified in microbiomes, it would break
272
this familiar conceptual connection between causation and corre
lation. Consider the region of
273
hysteresis (points S1-S3 and S5-S7, Fig. 4). The observed stead
y-state levels of
Kp
no longer
274
correlate with the input glucose
concentration. At 2 mM input
glucose, the system could be in
275
either the
Kp
-only state S3 or the
Kp–Bt
state S5. At ~650x10
6
CFU/mL of
Kp
, the input glucose
276
levels could be either 0.25 mM or 2 mM. Although there is no correlation between species’
277
abundance and input glucose concen
tration (in other words, know
ing glucose concentration is not
278
sufficient for predicting species
abundance; instead, a system’
s history must also be known), input
279
glucose concentration remains the causal factor. Furthermore,
under MSH, establishing causation
280
is insufficient for achieving
control: although input glucose c
oncentration is the causal factor
281
responsible for changes in the community state, it cannot be us
ed to fully control the community
282
(i.e. one cannot use changes in
glucose inputs to revert the
Kp–Bt
state back to the
Kp
-only state).
283
Alternative control strategies
(e.g. changes in oxygen levels o
r disruption of metabolic coupling),
284
derived from appropriate models, would need to be deployed unde
r MSH. Therefore, recognizing
285
whether and when MSH exists in human microbiomes will be critic
al for interpreting correlation
286
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
11
of
20
and causation, and for designing th
erapeutic control strategies
that can steer microbial communities
287
to desirable states.
288
289
290
291
292
Materials and Methods
293
Model development
294
For the computational simulations, we used the dynamic multispe
cies metabolic modeling
295
(DMMM) framework (
38
), which is an extension of dynamic flux balance analysis (dFBA
) applied
296
to microbial communities. In brief, the DMMM framework has two
components:
(Component 1)
297
external differential equations t
hat describe mass balances for
species and metabolite
298
concentrations in the CSTR (shown in Fig 2A and described in th
e SI). Unlike the traditional
299
method for solving differential equations in a bacterial system
, we do not assume that parameters
300
such as growth rates and metabolic flux rates are constant. Ins
tead, we allow the parameters to be
301
dynamic because we are studying a system with potentially rich
dynamics. In order to find the
302
values for these dynamic parameters we use flux balance analysi
s (
Component 2
) to solve for the
303
parameter values at every time step of the simulated time perio
d. Component 2 includes the
304
genome-scale models for each species. These models are used to
perform flux balance analysis at
305
every time point to obtain update
d parameters for the different
ial equations in Component 1.
306
Component 1:
The system is described as a continuous stirred tank reactor (C
STR) with the
307
following mathematical formulation:
308
ௗ
ௗ௧
ൌ휇
푋
െ
ி
ೠ
(1)
309
310
ௗௌ
ೕ
ௗ௧
ൌ
∑
푣
푋
ி
ௌ
ೕ
ିி
ೠ
ௌ
ೕ
(2)
311
312
ௗௌ
ೣ
ௗ௧
ൌ
∑
푣
௫௬
푋
퐾
푎ሺ푆
∗
െ푆
௫௬
ሻ
(3)
313
314
Here, V, is the volume of the reactor (constant),
푋
is the biomass (g/L) of the i
th
microbial species.
315
푆
is the concentration (mM) of the j
th
metabolite,
퐹
is the rate of flow (L/h) into the reactor,
퐹
௨௧
316
is the rate of flow (L/h) out of the reactor,
푆
ௗ
is the concentration of the j
th
metabolite in the feed
317
stream,
휇
(h
-1
) is the growth rate of the i
th
microbial species,
푣
is the metabolic flux of the j
th
318
substrate in the i
th
microbial species,
퐾
푎
is the volumetric oxygen transfer coefficient, and S* is
319
the dissolved oxygen satur
ation concentration.
320
In continuous culture, substrate
utilization can deviate from d
iauxic growth (as typically observed
321
in batch culture) and co-substrate utilization is possible
(
45
)
. The set of differential equations are
322
solved using the following analytical approximation:
323
324
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
12
of
20
푋
ൌ푋
,
푒
ቀఓ
ି
ಷ
ೇ
బ
ቁ∆்
(4)
325
326
푆
ൌ푆
∑
ሾ푣
బ
ఓ
బ
ିி
ሺ푋
,
푒
ቀఓ
ି
ಷ
ೇ
బ
ቁ∆்
െ푋
,
ሻሿ
ி
ሺௌ
ೕ
ିௌ
ೕ
ሻ
బ
∆푇
(5)
327
328
푆
௫௬
ൌ푆
௫௬
∑
ሾ푣
௫௬
బ
ఓ
బ
ିி
ሺ푋
,
푒
ቀఓ
ି
ಷ
ೇ
బ
ቁ∆்
െ푋
,
ሻሿ퐾
푎ሺ푆
∗
െ푆
௫௬
ሻ∆푇
(6)
329
330
At the beginning of every time step (
∆푇ሻ
, the parameters
휇
and
푣
are calculated using flux balance
331
analysis (FBA) from genome-scale models (Component 2) and fed back into equations (4), (5), and
332
(6). This process is repeated for all time intervals in the sim
ulated time period.
333
334
Component 2:
Genome-scale metabolic models ar
e used to establish genotype-p
henotype
335
relationships and capture the metabolic capabilities of each mo
del organism. Furthermore, these
336
models allow us to integrate me
tabolic network topology information with RNA sequencing data
337
for transcriptomic analysis of the CSTR experiments (described
further in “RNA Sequencing and
338
Analysis” section). We used the published iYL1228 model of
Klebsiella pneumoniae (Kp)
MGH
339
78578 (
37
) and the published
iAH991 model of
Bacteroides thetaiotaomicron (Bt)
VPI 5482
(
36
)
.
340
Because these models are well validated (and curated), we added
the minimum number of
341
parameters that allow for integration of these models to the co
mmunity dynamic flux balance
342
analysis framework. Our changes
include adding a pathway for de
xtran uptake and hydrolysis to
343
glucose in the
Bt
iAH991 model. The pathway lumps hydrolysis of dextran to glucose into a single
344
reaction. In this lumped reaction we assume that 50% of the glucose produced from dextran by
Bt
345
can be released into the environment for shared use. For the purpose of the simulations, dextran is
346
assumed to be 100 glucose units. The genome-scale models are so
lved separately for each species
347
by flux balance analysis (FBA) (
46
) at each time point:
348
349
max푐
்
푣
ప
ഥ
(7)
350
푠.푡. 퐴
푣
ప
ഥൌ0
351
푣̅
,
൏푣
ప
ഥ൏ 푣 ̅
,௨
352
353
Where c is the cost vector,
푣̅
is the vector of fluxes, and
A
is the matrix of mass balance
354
stoichiometries. The optimization criterion is biomass growth r
ate (for each species). For the
355
bounds for the fluxes, we used the values in the curated, publi
shed genome-scale models (
Kp
356
iYL1228 and
Bt
iAH99). The uptake fluxes explicitly modeled in Component 1 (d
extran, glucose,
357
acetate, oxygen) are bounded by Mic
haelis–Menten kinetics:
358
푣
,௨
ൌ푣
,௫
ௌ
ೕ
ାௌ
ೕ
(8)
359
The values for
푣
,௫
and
퐾
for some of the metabolites in the model were estimated from b
atch
360
experiments. Batch culture experiments were carried out in a 96
-well flat-bottom plate. Overnight
361
cultures grown anaerobically in minimal medium supplemented wit
h either 0.5% w/v dextran or
362
0.5% w/v glucose were diluted 1:20 (for
Bt
) and 1:100 (for
Kp
), and outgrown to mid-log phase.
363
The cultures were then pelleted and re-suspended at OD 1 (for
Bt
) and OD 0.1 (for
Kp
) in carbon-
364
free minimal medium. We added 10 μL of cells to 200 μL of minim
al medium containing various
365
concentrations (0.125 – 0.5% w/v) of the carbon source. The pla
te was incubated at 37 °C and
366
OD600 was measured every 10 min. For batch cultures, Monod growth kinetics was assumed:
367
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
13
of
20
ௗ
ௗ௧
ൌ푋휇
௫
ௌ
ௌା
(9)
368
ௗௌ
ௗ௧
ൌെ푣
௫
ௌ
ௌା
(10)
369
Where
휇
௫
is the maximum growth rate the given bacteria can achieve on a
carbon source when
370
it is not resource limited. Growth data from replicate wells of
multiple concentrations of carbon
371
source were fitted simultaneously using Bayesian parameter esti
mation implemented with Markov
372
chain Monte Carlo (MCMC) (
47
). Individual growth curves were allowed to have distinct initi
al
373
cell concentrations and background values, with other parameters held constant. The fitted
374
parameters are presented in Table S1 and Fig S3. For all other
metabolites captured in the
375
differential equations, the K
m
and v
max
values are assumed to be the same (v
max
of 10
376
mmol/gCDW
h and K
m
of 0.01 mM), based on literature for
Escherichia coli
(
48, 49
)
.
377
Values for parameters and initial conditions used in the model
are presented in Table S2. Initial
378
conditions are chosen to represen
t the experimental setup, wher
eby we first establish a steady state
379
for
Kp
in the CSTR before inoculating
Bt
. Therefore, in the models we start with a higher
380
concentration of
Kp
than
Bt
. The initial conditions for
Kp
in the reactor is arbitrarily chosen to be
381
the experimentally measured mono-c
ulture steady concentration o
f
Kp
at an input glucose
382
concentration of 0.25 mM. Th
e initial conditions for
Bt
in the reactor is 0.0015 g/L, which is
383
equivalent to addition of 1 mL of OD 1
Bt
into the reactor, as done expe
rimentally. For most steady
384
state conditions, glucose is lim
iting and therefore the initial
conditions for glucose concentration in
385
the reactor is chosen to be 0mM.
386
We note that although an individua
l FBA simulation (Component 1
) is
387
linear, a dynamic FBA (the combin
ation of Component 1 and Compo
nent 2) is not linear. The
388
behavior of the cells will shift dramatically (a nonlinear resp
onse) when substrate levels in the
389
environment cross certain thresholds. The two biggest sources o
f nonlinearity in our model
390
system are binary growth/no-growth behavior of
Bt
in the presence or absence of O
2
, and the binary
391
capacity of
Kp
to utilize carbon from dextran (through its breakdown into oligosaccharides by
Bt
)
392
in the presence or absence of
Bt
. To model the binary growth/no-growth behavior of
Bt
in the
393
presence of O
2
, we included a condi
tional operator before the FBA simulations
for each time step:
394
If the O
2
levels are above 350 nM,
the growth rate of
Bt
is set to zero and the FBA is not run for
Bt
.
395
Whereas if O
2
levels are lower than 350 nM O
2
, then the FBA simulation for
Bt
is allowed
396
to proceed. The value for the O
2
growth threshold for the model is arbitrarily chosen from the
397
concentration range repor
ted for the aerotolerant
Bt
(
50
).
398
To computationally identify the regions of stability with respe
ct to glucose and oxygen (Fig. 3C
399
and Fig. S4), we varied oxygen i
nput flow rates at constant inp
ut glucose concentration for each
400
glucose condition examined. We evaluated 11 glucose conditions
ranging from 1 mM to 6 mM. For
401
each given glucose input concentration, we started with oxygen
at an input flow rate of 6 mL/min
402
and ran the simulation for 50 h to ensure the system reached a
steady state. We then decreased the
403
oxygen input by 0.5 mL/min intervals down to 0.5 mL/min; for ea
ch oxygen condition, ensuring
404
the system reached a steady state (we refer to the 6mL/min - 0.
5mL/min oxygen variations for a
405
given constant glucose input con
centration as the “forward simu
lations”). The concentration of
406
oxygen input at which
Bt
starts to grow is identified as the “tipping point” to the mon
o-stable
Kp–
407
Bt
state. After running the 0.5mL/min oxygen simulation we increa
sed the concentration back to 6
408
mL/min at intervals of 0.5 mL/min (we refer to the 0.5 mL/min - 6 mL/min oxygen variations as
409
the “reverse simulations”). The
concentration of oxygen at whic
h
Bt
can no longer grow and gets
410
washed out is identified as “tipping point” to the mono-stable
Kp-
only state. The region between
411
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
14
of
20
these two tipping points to the mono-stable
Kp–Bt
state in the forward simulations and the mono-
412
stable
Kp
-only state in the reverse simulations is identified as the reg
ion of bi-stability. The colors
413
in Fig. S4 represent the steady state concentration of
Kp
in the “reverse simulations” divided by the
414
steady state concentration of
Kp
in the “forward simulations.” In regions of mono-stability the
415
concentration of
Kp
is the similar in both the “forward” and “reverse simulations”
, and therefore
416
has a value of approximately 1.
417
Continuous culture of
K. pneumoniae
and
B. thetaiotaomicron
418
Continuous culture experiments w
ere carried out in a 500 mL bio
reactor (Mini-bio Applikon
419
Biotechnology, Delft, Netherlands) with a total culture volume
of 200 mL. Minimal media (3.85
420
g/L KH
2
PO
4
, 12.48 g/L K
2
HPO
4
, 1.125 g/L (NH
4
)
2
SO
4
, 1X MMS (20X MMS:
17.6 g/L NaCl, 0.4
421
g/L CaCl
2
, 0.4 g/L MgCl
2
×6H2O, 0.2 g/L MnCl
2
×4H2O, 0.2 g/L CoCl
2
×6H2O), 10 mL/L Wolfe’s
422
mineral solution (
51
), 10 mL/L Wolfe’s vitamin solution (
51
), 4.17 μM FeSO4×7H20, 0.25 mM
423
cysteine, 1 μM menadione, 2 μM r
esazurin, 1 g/L
dextran (Sigma
D5376, avg. mol. wt 1.5e6-2e6)
424
and glucose at varying concentrations) was purged with 100% N
2,
stored under anaerobic conditions
425
prior to use, and maintained under N
2
during operation of the CSTR. We calibrated the dissolved
426
oxygen probe by aerating the reactor with CO
2
(5 ml/min) and N
2
(45 ml/min). The stable
427
measurement without O
2
input was taken to be 0 μM dissolved oxygen. We then sparged t
he reactor
428
with O
2
(1.7 ml/min), CO
2
(5ml/min), and N
2
(43.3 ml/min). This stable measurement was taken to
429
be 30.714 μM O
2
(calculated assuming 6.056 mg/L
dissolved oxygen at 1 atm and 37 °C).
430
431
For the experiments, the bioreactor was sparged with 50 mL/min
total gas (1.7 mL/min O
2
, 5
432
mL/min CO
2
, and balance of N
2
), and agitated with two six-bladed Rushton turbines operated a
t
433
750 rpm. Temperature was maintai
ned at 37°C, and
a residence time of 5 h (input and output flow
434
rates of 40 ml/h flowrate) was used for all experiments. Dissol
ved oxygen, pH, and biomass were
435
monitored throughout. For in
itial inoculation of
Kp
, 1 mL of OD 1 culture was injected through the
436
septum, and grown in batch cultu
re until stationary phase (indi
cated by a levelled biomass readout
437
and an increase of dissolved oxygen levels) before beginning co
ntinuous culture. For every
438
experimental condition examined (input glucose concentration) i
n the CSTR, we waited for the
439
system to first reach steady state (at least 24 h). To assess w
hether a steady state had been reached,
440
we monitored the total biomass in the reactor using a real-time
OD probe. Once the system reached
441
steady state, we took three samples over the course of 24-48 h. The time interval between each
442
sample collection was at least
one residence time (5 h). Residence time is defined as the time
it
443
takes to entirely exchange the volume of the reactor.
For introduction of
Bt
, a log phase (OD 0.6-
444
0.8) anaerobic culture grown in minimal media with 0.5% dextran
and 2 mM cysteine was pelleted
445
(5 min at 3500 g) and washed twice using dextran/glucose-free anaerobic minimal media. Cells
446
were carbon-starved at 37 °C fo
r 30 min, washed (once), and re-
suspended in dextran/glucose-free
447
minimal media to OD 1. We used 1 mL of this
Bt
cell suspension for inocul
ation into the reactor
448
and a sample was collected immediately after inoculation. A sub
sequent sample was collected for
449
quantification after at least 2 residence times had passed. In
the
Kp
-only state conditions (0.25 mM,
450
1 mM, 2 mM glucose),
Bt
is washed out, as described in the results section. To ensure
451
reproducibility of a washout for these conditions, the
Bt
inoculation and sample collection process
452
was repeated a total of three times. In
Kp–Bt
state conditions (5 mM, 2
mM, 1 mM, 0.25 mM, and
453
0 mM glucose), where
Bt
growth persisted, re-inoculation of
Bt
was no longer necessary for each
454
new glucose steady state condition; three samples separated by
at least one residence time were
455
collected for each steady state condition. To collect samples, ~0.5 mL of culture was removed from
456
the bioreactor in a 3 mL luer-lock syringe and discarded before collection of 1.5–2 mL culture.
457
Supernatant from 700 μL of the co
llected sample was stored at -80 °C for SCFA analysis, a 50 μL
458
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
15
of
20
sub-sample was treated with DNAse (2.5 μL of NEB DNase I 2000 u
/mL per 50 μL) for subsequent
459
DNA extraction, and two 250 μL aliquot
s were used for extractio
n of RNA.
460
461
Quantification of bacterial abundance
462
463
CSTR culture samples were treat
ed with NEB DNAse I (100 u/mL fi
nal concentration) for 10 min
464
at 37 °C immediately after collection. DNA was extracted using
the ZyGEM prepGEM Bacteria kit
465
(ZyGEM, Southampton, England) according to the manufacturer pro
tocol. Samples were extracted
466
in 100 uL total volume (20 μL culture sample and 80 μL of extraction mixture), incubated at 37 °C
467
for 15 min, 75 °C for 5 min, 95 °C for 5 min, then cooled to 4
°C. DNA was stabilized by adding
468
10X TE to a final concentration of
1X TE before storage at 4 °C.
469
470
qPCR quantification
Extracted DNA was quantified by qPCR using the Eco Real-time PC
R
471
system (Illumina, San Diego, CA, USA). The components in the qP
CR mix used in this study were
472
as follows: 1 μL of extracted DNA, 1X SsoFast™ EvaGreen Supermi
x (Bio-Rad Laboratories,
473
Hercules, CA, USA), 500 nM forward primer, and 500 nM reverse p
rimer. For detection of each
474
bacterial species in the community primer sets specific to
Bt
(forward primer: 5′-
475
GGAGTTTTACTTTGAATGGAC-3′; rever
se primer: 5′-CTGCCCTTTTACAATGGG-3’) and
476
Kp
(forward primer: 5′-ATTTGAAGAGGTTGCAAACGAT-3′; reverse primer:
5′-
477
TTCACTCTGAAGTTTTCTTGTGTT-3′) were
used. Quantification of cell
concentrations were
478
determined using DNA standards of single species prepared using 10X serial dilutions of log phase
479
cultures extracted as above. Cell concentrations of standards w
ere determined by hemocytometer.
480
For conversion of OD and cell c
oncentration to biomass concentr
ation (gram cell dry weight/L),
481
100 mL of culture for each individual species incubated anaerob
ically at 37 °C was harvested and
482
pellets were dried at 80°C for
~48 h before recording mass.
483
484
Digital PCR quantification
Archived DNA samples from the CSTR were quantified by digital
485
PCR (dPCR) using a QX200 Droplet Digital PCR System (Bio-Rad).
The components in the dPCR
486
mix were as follows: 1 μL of dilutions of extracted DNA, 1X QX200 ddPCR EvaGreen Supermix
487
(Bio-Rad), 500 nM forward primer, and 500 nM reverse primer. Fo
r detection of each bacterial
488
species in the community primer sets specific to
Bt
(forward primer:5′-
489
GGTGTCGGCTTAAGTGCCAT-3′; reverse
primer: 5′-CGGAYGTAAGGGCCGTGC-
3′) and
Kp
490
(forward primer:5′-ATGGCTGTCGTCAGCTCGT-3′; reverse primer: 5′-
491
CCTACTTCTTTTGCAACCCACTC-3′) were used. The dPCR mix was loaded
into DG8
492
Cartridges (Bio-Rad), which were filled with QX200 Droplet Generation Oil for EvaGreen (Bio-
493
Rad) and loaded into the QX200 Droplet Generator (Bio-Rad). The
generated droplets were
494
transferred to the C1000 Touch Thermal Cycler (Bio-Rad) for the
following protocol: 5 min at 95
495
°C, 40 cycles of 30 sec at 95 °C, 30 sec at 60 °C (
Kp
) or 65 °C (
Bt
), and 30 sec at 72 °C (ramping
496
rate reduced to 2°C/s), and final dye stabilization steps of 5
min at 4 °C, and 5 min at 90 °C. The
497
stabilized plates were loaded into the QX200 Droplet Reader and
analyzed using the QuantaSoft
498
Analysis Software (Bio-Rad).
499
500
RNA sequencing and analysis
501
From the CSTR samples, a 250 μL aliquot was used for metatransc
riptomic analysis. The freshly
502
collected CSTR sample was immedi
ately placed into Qiagen RNAprotect Bacteria Reagent
503
(Qiagen, Hilden, Germany) for RNA stabilization. RNA was extracted using the Enzymatic Lysis
504
of Bacteria protocol of the Qia
gen RNeasy Mini Kit and processed according to the manufacturer’s
505
protocol. DNA digestion was perf
ormed during extraction using t
he Qiagen RNase-Free DNase
506
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
16
of
20
Set. The quality of extracted RNA was measured using an Agilent
2200 TapeStation (Agilent, Santa
507
Clara, CA, USA). Extracted RNA samples were prepared for sequen
cing using the NEBNext Ultra
508
RNA Library Prep Kit for Illumina (New England Biolabs, Ipswitch, MA, USA) and the
509
NEBNeE
xt Multiplex Oligos for Illumina. Libraries were sequenced at 1
00 single base pair reads
510
and a sequencing depth of 10 million reads on an Illumina HiSeq
2500 System (Illumina, San
511
Diego, CA, USA) at the
Millard and Muriel J
acobs Genetics and G
enomics Laboratory, California
512
Institute of Technology. Raw reads from the sequenced libraries
were subjected to quality control
513
to filter out low-quality reads and trim the adaptor sequences
using Trimmomatic (v. 0.35). Because
514
our samples were a mixture of
Kp
and
Bt
cells, to separate the reads for each species, we did the
515
following: Reads that aligned to rRNA and tRNA of
Bt
and
Kp
were first removed, as those
516
sequences contain overlapping reads between the two species. Each sample was then separately
517
aligned to
Bt
VPI-5482 (Genome accession number: GCA_000011065.1) and
Kp
MGH-78578
518
(Genome accession number: GCA_000016305.1) using Bowtie2 (v. 2.
2.5) and quantified using the
519
Subread package (v. 1.5.0-p1). Ge
ne expression was defined in t
ranscripts per million (TPM) for
520
each species and gene expression analysis was performed using D
ESeq2 (v. 1.22.2; default settings,
521
which provides two-ta
ilored p-values).
522
523
To determine the most differentia
lly regulated metabolic pathwa
y between the
Kp-Bt
and
Kp
-only
524
states we used an approach by Patil et al. (
41
), which combines gene expression data with
525
topological information collected from genome-scale metabolic m
odels. In brief, in this approach
526
the metabolic network is presented as a bipartite undirected gr
aph, where metabolites as well as
527
enzymes are represented as nodes
(this graph is obtained from g
enome-scale models). Differential
528
data can be mapped on the enzyme nodes of the graph with specif
ication of the significance of
529
differential gene expression for each enzyme, i. We used DESeq
to perform our differential gene
530
expression analysis between sample points in the two states (
Kp
-only and
Kp-Bt
) to obtain
P
-values.
531
The
P
-values are subsequently converte
d to Z-score for an enzyme nod
e using the inverse normal
532
cumulative distribution. Finally, the Z-score of each metabolit
e node is calculated based on the
533
normalized transcriptional res
ponse of its k neighboring enzyme
s:
534
535
푍
௧௧
ൌ
1
√
푘
푍
ଵ
536
The metabolites with the highest Z-scores mark the pathways wit
h substantial regulation between
537
two states.
538
539
540
541
542
543
544
Supplementary Materials
545
546
Fig. S1.
The workflow for the continuous bioreactor experiment
547
Fig. S2.
Imaging of samples collected fro
m the continuously stirred tan
k reactor experiments.
548
Fig. S3.
Bayesian parameter fitting of K and v
max
to experimental batch growth data.
549
Fig. S4.
A quantitative view of regions of stability as a function of glucose concentrations in the
550
input feed and oxygen flow rate
s into the reactor.
551
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint
Page
17
of
20
Fig S5.
State switch from
Kp
-only state (state A) to the
Kp–Bt
state (state B) by the transient
552
addition (pulse) of shor
t chain fatty acids.
553
Fig S6.
Quantification of
B. thetaiotaomicron
and
K. pneumoniae
in archived CSTR samples using
554
digital PCR
555
Fig. S7.
Simulations of extracellular c
oncentrations of dextran and glu
cose (A), and acetate and
556
oxygen (B) as a factor of input g
lucose concentration under con
stant oxygen feed (
flow rate of 1.7
557
mL/min).
558
Fig. S8.
P
-values and individual data point
s for (A) Fig. 3A (B) Fig. 3B
(C) Fig. S5 and (D) Fig.
559
S6 of the manuscript.
560
Fig. S9.
A pilot CSTR experiment demons
trating state-switching and hysteresis.
561
Table S1.
Bayesian parameter estimation for K and v
max
used in the Michaelis–Menten equations
562
to constrain nutrient uptake fl
ux rates for flux balance analys
is calculations.
563
Table S2.
Values for the parameters use
d in the dynamic flux balance ana
lysis simulations.
564
Table S3.
The top scoring 50 metabolites involved in the most regulated
metabolic pathways,
565
ordered by Z-score value.
566
Table S4.
P
-values between statistically significant samples (gene express
ion values) for K.
567
pneumoniae in Fig. 5.
568
569
570
571
References and Notes
572
573
1. V. Singh, S. Proctor, B. Willing, Koch's postulates, microbi
al dysbiosis and inflammatory
574
bowel disease.
Clin. Microbiol. and Infect.
22
, 594-599 (2016).
575
2. L. Zhao, The gut microbiota and obesity: from correlation to
causality.
Nat. Rev.
576
Microbiol.
11
, 639 (2013).
577
3. R. Blumberg, F. Powrie, Micro
biota, disease, and back to health: a metastable journey.
578
Sci. Transl. Med.
4
, 137rv137 (2012).
579
4. J. Bures
et al.
, Small intestinal bacterial overgrowth syndrome.
World J. Gastroenterol.
580
16
, 2978-2990 (2010).
581
5. L. Zheng, C. J. Kelly, S. P. C
olgan, Physiologic hypoxia and
oxygen homeostasis in the
582
healthy intestine. A Review in t
he Theme: Cellular Responses to Hypoxia.
American
583
journal of physiology. Cell physiology
309
, C350-C360 (2015).
584
6. D. Bradshaw, P. Marsh, G. Watson, C. Allison, Oral anaerobes
cannot survive oxygen
585
stress without interacting with f
acuItative/aerobic species as
a microbial commmunity.
586
Lett. Appl. Microbiol.
25
, 385-387 (1997).
587
7. D. Gonze, L. Lahti, J. Raes, K
. Faust, Multi-stability and the origin of microbial
588
community types.
ISME J.
11
, 2159 (2017).
589
8. S. A. Shetty, F. Hugenholtz, L. Lahti, H. Smidt, W. M. de Vo
s, Intestinal microbiome
590
landscaping: insight in commun
ity assemblage and implications f
or microbial modulation
591
strategies.
FEMS Microbiol. Rev.
41
, 182-199 (2017).
592
9. V. Bucci, S. Bradde, G. Biroli
, J. B. Xavier, Social interaction, noise and antibiotic-
593
mediated switches in the
intestinal microbiota.
PLoS Comp. Biol.
8
, e1002497 (2012).
594
10. S. Vet
et al.
, Bistability in a system of tw
o species interacting through mu
tualism as well
595
as competition: Chemostat vs. Lotka-Volterra equations.
PloS one
13
, e0197462 (2018).
596
11. D. R. Amor, C. Ratzke, J. Gor
e, Transient invaders can indu
ce shifts between alternative
597
stable states of microbial communities.
Science Advances
6
, eaay8676 (2020).
598
12. V. L. Andersen Woltz, C. I. Abreu, J. Friedman, J. Gore, Emergence of alternative stable
599
states in microbial communities in a fluctuating environment.
bioRxiv
, 678367 (2019).
600
.
CC-BY 4.0 International license
author/funder. It is made available under a
The copyright holder for this preprint (which was not peer-reviewed) is the
.
https://doi.org/10.1101/2020.02.28.968941
doi:
bioRxiv preprint