Processes 2019, 7,97
data [12,61–64] and have reduced the glucose concentration and increased AA concentrations for the
dysbiosis case consistent with experimental observation [33,110]. We performed a sensitivity analysis
of these concentrations to show that a similar behavior (i.e., healthy state) as that reported for the
nominal values occurred if the glucose to AA ratio was sufficiently large (Figure S2). By contrast, a
CDI dysbiosis-like state was obtained when the glucose to AA ratio was sufficiently small.
Uptake rates of nutrients and byproducts were assumed to follow Michaelis–Menten kinetics.
Due to lack of available data, maximum uptake rates and Michaelis–Menten constants were assumed
to be independent of species and metabolite. Calculated uptake rates were imposed as lower bounds of
the exchange fluxes in the species metabolic reconstructions. The calculated growth rate, uptake fluxes,
and secretion fluxes from each reconstruction served as inputs to reaction-diffusion-type equations for
the biomass concentration of each species and the molar concentration of each nutrient and byproduct.
This formulation yielded a set of 23 partial differential equations (PDEs) in the time and the axial
direction z with embedded linear programs (LPs) for species metabolism (see Appendix S1). Following
our previous methodology [50,51], lexicographic optimization with growth rate maximization as
the primary objective was used to avoid alternative optima that would render the biofilm model
non-smooth. This approach yielded a total of 71 LPs.
The biofilm model equations were solved by spatially discretizing the PDEs into a large set
of ordinary differential equations (ODEs) [111,112]. We used 25 spatial node points to achieve a
suitable compromise between solution accuracy and computational efficiency, which produced a
discretized model with 575 ODEs and 1775 LPs that was solved with the MATLAB code DFBAlab [113].
We used Gurobi 6.5.2 for the LP solution, the stiff MATLAB solver ode15s for ODE integration,
and DFBAlab running in MATLAB 9.0 (R2016a). Although not explored here, our biofilm modeling
method can be extended to more species and extracellular metabolites. For N spatial discretization
points, the addition of each new extracellular metabolite would generate N additional ODEs. For m
total extracellular metabolites, the addition of each new species would generates N additional ODEs
and m + 1 LPs. Because the LP solution scales more favorably than the ODE solution, we anticipated
that models with approximately 1000 ODEs and 7500 LPs would remain computationally viable on a
typical desktop computer. These equation numbers translate into approximately 10 species and 30
extracellular metabolites.
4.2. Biofilm Model Parameterization and Tuning
Nominal parameter values used in the multispecies biofilm model are shown in Table 2.
The parameters were obtained from the experimental literature to the extent possible and from
our previous modeling studies [50,51] as necessary. The bulk glucose and amino acid concentrations
at the biofilm–stool interface were specified to reflect healthy gut conditions. Due to the lack of
species-specific uptake data, we used published kinetic parameters reported for E. coli [114]. Due to
the lack of data, all eight byproducts were assumed to have the same uptake parameters as glucose.
For simplicity, all eight amino acids were assumed to have the same uptake parameters obtained as
the average of amino acid-dependent values reported for E. coli [114].
With all other parameter values fixed, the biofilm model was qualitatively tuned to achieve
biomass and SCFA fractions within experimental ranges for a healthy patient. The species abundances
were tuned by adjusting the non-growth-associated ATP maintenance (ATPM) values of the four
metabolic reconstructions following our previous studies [50,51]. Our justification for tuning these
values was the simple nature of the biofilm model, which neglected other phyla (e.g., Actinobacteria),
other nutrients (e.g., oligosaccharides, fats), other species interactions (e.g., Actinobacteria cross-feeding
of SCFAs and organic acids), as well as host metabolism present in the actual gut environment.
These ATPM values listed in Table 2 produced B. thetaiotaomicron:F. prausnitzii:E. coli:C. difficile
abundances of 71%:21%:7%:1%, which were deemed reasonable based on published data [57,58].
We found that the coexistence of the four species was achieved over a range of ATPM values
(not shown here).
32
data [12,61–64] and have reduced the glucose concentration and increased AA concentrations for the
dysbiosis case consistent with experimental observation [33,110]. We performed a sensitivity analysis
of these concentrations to show that a similar behavior (i.e., healthy state) as that reported for the
nominal values occurred if the glucose to AA ratio was sufficiently large (Figure S2). By contrast, a
CDI dysbiosis-like state was obtained when the glucose to AA ratio was sufficiently small.
Uptake rates of nutrients and byproducts were assumed to follow Michaelis–Menten kinetics.
Due to lack of available data, maximum uptake rates and Michaelis–Menten constants were assumed
to be independent of species and metabolite. Calculated uptake rates were imposed as lower bounds of
the exchange fluxes in the species metabolic reconstructions. The calculated growth rate, uptake fluxes,
and secretion fluxes from each reconstruction served as inputs to reaction-diffusion-type equations for
the biomass concentration of each species and the molar concentration of each nutrient and byproduct.
This formulation yielded a set of 23 partial differential equations (PDEs) in the time and the axial
direction z with embedded linear programs (LPs) for species metabolism (see Appendix S1). Following
our previous methodology [50,51], lexicographic optimization with growth rate maximization as
the primary objective was used to avoid alternative optima that would render the biofilm model
non-smooth. This approach yielded a total of 71 LPs.
The biofilm model equations were solved by spatially discretizing the PDEs into a large set
of ordinary differential equations (ODEs) [111,112]. We used 25 spatial node points to achieve a
suitable compromise between solution accuracy and computational efficiency, which produced a
discretized model with 575 ODEs and 1775 LPs that was solved with the MATLAB code DFBAlab [113].
We used Gurobi 6.5.2 for the LP solution, the stiff MATLAB solver ode15s for ODE integration,
and DFBAlab running in MATLAB 9.0 (R2016a). Although not explored here, our biofilm modeling
method can be extended to more species and extracellular metabolites. For N spatial discretization
points, the addition of each new extracellular metabolite would generate N additional ODEs. For m
total extracellular metabolites, the addition of each new species would generates N additional ODEs
and m + 1 LPs. Because the LP solution scales more favorably than the ODE solution, we anticipated
that models with approximately 1000 ODEs and 7500 LPs would remain computationally viable on a
typical desktop computer. These equation numbers translate into approximately 10 species and 30
extracellular metabolites.
4.2. Biofilm Model Parameterization and Tuning
Nominal parameter values used in the multispecies biofilm model are shown in Table 2.
The parameters were obtained from the experimental literature to the extent possible and from
our previous modeling studies [50,51] as necessary. The bulk glucose and amino acid concentrations
at the biofilm–stool interface were specified to reflect healthy gut conditions. Due to the lack of
species-specific uptake data, we used published kinetic parameters reported for E. coli [114]. Due to
the lack of data, all eight byproducts were assumed to have the same uptake parameters as glucose.
For simplicity, all eight amino acids were assumed to have the same uptake parameters obtained as
the average of amino acid-dependent values reported for E. coli [114].
With all other parameter values fixed, the biofilm model was qualitatively tuned to achieve
biomass and SCFA fractions within experimental ranges for a healthy patient. The species abundances
were tuned by adjusting the non-growth-associated ATP maintenance (ATPM) values of the four
metabolic reconstructions following our previous studies [50,51]. Our justification for tuning these
values was the simple nature of the biofilm model, which neglected other phyla (e.g., Actinobacteria),
other nutrients (e.g., oligosaccharides, fats), other species interactions (e.g., Actinobacteria cross-feeding
of SCFAs and organic acids), as well as host metabolism present in the actual gut environment.
These ATPM values listed in Table 2 produced B. thetaiotaomicron:F. prausnitzii:E. coli:C. difficile
abundances of 71%:21%:7%:1%, which were deemed reasonable based on published data [57,58].
We found that the coexistence of the four species was achieved over a range of ATPM values
(not shown here).
32