# Hydrogenotrophic methanogens of the mammalian gut: Functionally similar, thermodynamically different-A modelling approach

Rafael Munoz Tamayo, Milka Popova, Maxence Tillier, Diego Morgavi, Jean-Pierre Morel, Gérard Fonty, Nicole Morel-Desrosiers

<figure><img src="image_1.png" /></figure>

## To cite this version:

Rafael Munoz Tamayo, Milka Popova, Maxence Tillier, Diego Morgavi, Jean-Pierre Morel, et al. Hydrogenotrophic methanogens of the mammalian gut: Functionally similar, thermodynamically different-A modelling approach. PLoS ONE, 2019, 14 (12), pp.e0226243. 10.1371/journal.pone.0226243. hal-02473686

HAL Id: hal-02473686

https://hal.science/hal-02473686

Submitted on 10 Feb 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Hydrogenotrophic methanogens of the mammalian gut: Functionally similar, thermodynamically different—A modelling approach

Rafael Muñoz-Tamayo<sup>1,*</sup>, Milka Popova<sup>2,†</sup>, Maxence Tillier<sup>2</sup>, Diego P. Morgavi<sup>2,†</sup>, Jean-Pierre Morel<sup>3</sup>, Gérard Fonty<sup>3</sup>, Nicole Morel-Desrosiers<sup>3</sup>

1 UMR Modélisation Systémique Appliquée aux Ruminants, INRA, AgroParisTech, Université Paris-Saclay, Paris, France, 2 Institute National de la Recherche Agronomique, UMR1213 Herbivores, Clermont Université, VetAgro Sup, UMR Herbivores, Clermont-Ferrand, France, 3 Université Clermont Auvergne, CNRS, LMGE, Clermont-Ferrand, France

© These authors contributed equally to this work.

* rafael.munoz-tamayo@inra.fr

<figure><img src="image_3.png" /></figure>

## Abstract

Methanogenic archaea occupy a functionally important niche in the gut microbial ecosystem of mammals. Our purpose was to quantitatively characterize the dynamics of methanogenesis by integrating microbiology, thermodynamics and mathematical modelling. For that, *in vitro* growth experiments were performed with pure cultures of key methanogens from the human and ruminant gut, namely *Methanobrevibacter smithii*, *Methanobrevibacter ruminantium* and *Methanobacterium formicium*. Microcalorimetric experiments were performed to quantify the methanogenesis heat flux. We constructed an energetic-based mathematical model of methanogenesis. Our model captured efficiently the dynamics of methanogenesis with average concordance correlation coefficients of 0.95 for CO₂, 0.98 for H₂ and 0.97 for CH₄. Together, experimental data and model enabled us to quantify metabolism kinetics and energetic patterns that were specific and distinct for each species despite their use of analogous methane-producing pathways. Then, we tested *in silico* the interactions between these methanogens under an *in vivo* simulation scenario using a theoretical modelling exercise. *In silico* simulations suggest that the classical competitive exclusion principle is inapplicable to gut ecosystems and that kinetic information alone cannot explain gut ecological aspects such as microbial coexistence. We suggest that ecological models of gut ecosystems require the integration of microbial kinetics with nonlinear behaviours related to spatial and temporal variations taking place in mammalian guts. Our work provides novel information on the thermodynamics and dynamics of methanogens. This understanding will be useful to construct new gut models with enhanced prediction capabilities and could have practical applications for promoting gut health in mammals and mitigating ruminant methane emissions.

<figure><img src="image_4.png" /></figure>

### OPEN ACCESS

**Citation:** Muñoz-Tamayo R, Popova M, Tillier M, Morgavi DP, Morel J-P, Fonty G, et al. (2019) Hydrogenotrophic methanogens of the mammalian gut: Functionally similar, thermodynamically different—A modelling approach. PLoS ONE 14 (12): e0226243. https://doi.org/10.1371/journal.pone.0226243

**Editor:** James E. Wells, United States Department of Agriculture, Agricultural Research Service, UNITED STATES

**Received:** July 11, 2019

**Accepted:** November 24, 2019

**Published:** December 11, 2019

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0226243

**Copyright:** © 2019 Muñoz-Tamayo et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
**Data Availability Statement:** All relevant data are
within the manuscript and its Supporting
Information files.

**Funding:** This work received funding from Inra
PHASE department and the Inra MEM
metaprogramme to MP, RMT. The funder had no
role in study design, data collection and analysis,
decision to publish, or preparation of the
manuscript.

**Competing interests:** The authors have declared
that no competing interests exist.

**Introduction**

Methanogenic archaea inhabit the gastro-intestinal tract of mammals where they have established syntrophic interactions within the microbial community [1–3] playing a critical role in the energy balance of the host [4,5]. In the human gut microbiota, the implication of methanogens in host homeostasis or diseases is poorly studied, but of growing interest [6]. *Methanobrevibacter smithii* (accounting for 94% of the methanogen population) and *Methanosphaera stadtmanae* are specifically recognized by the human innate immune system and contribute to the activation of the adaptive immune response [7]. Decreased abundance of *M. smithii* was reported in inflammatory bowel disease patients [8], and it has been suggested that methanogens may contribute to obesity [9]. In the rumen, the methanogens community is more diverse though still dominated by *Methanobrevibacter* spp., followed by *Methanomicrobium* spp., *Methanobacterium* spp. [10] and *Methanomassilicoccus* spp [11]. However, the proportion of these taxa could vary largely, with *Methanomicrobium mobile* and *Methanobacterium formicium* being reported as major methanogens in grazing cattle [12]. Though methanogens in the rumen are essential for the optimal functioning of the ecosystem (by providing final electron acceptors), the methane they produce is emitted by the host animal, contributing to global greenhouse gas (GHG) emissions. In the gastrointestinal tract of mammals, major rumen methanogens [13] and the dominant human archaeon *M. smithii* [14], are hydrogenotrophic archaea without cytochrome (membrane-associated electron transfer proteins). Cytochrome-lacking methanogens exhibit lower growth yields than archaea with cytochromes [15]. However, this apparent energetic disadvantage has been counterbalanced by a greater adaptation to the environmental conditions prevailing in the gastrointestinal tract [16], and by the establishment of syntrophic interactions with feed fermenting microbes. This syntrophic cooperation centred on hydrogen allows the anaerobic reactions of substrate conversion to proceed close to the thermodynamic equilibrium [17,18] (that is with Gibbs free energy change close to zero).

To our knowledge, thermodynamic considerations on human gut metabolism have been poorly addressed in existing mathematical models [19–22], although the theoretical frameworks have been developed in other domains to calculate stoichiometric and energetic balances of microbial growth from the specification of the anabolic and catabolic reactions of microbial metabolism [23,24], and advances have been done to link thermodynamics to kinetics [25–27]. For the rumen, thermodynamic principles have been incorporated already into mathematical research frameworks because of their important role in feed fermentation. Thermodynamic studies have been performed to investigate theoretically (i) the profile of fermentation [28], (ii) alternative routes for hydrogen utilization [29], and (iii) the effect of hydrogen partial pressure on glucose fermentation and methanogenesis [30,31]. However, studies assessing quantitative comparisons between experimental data and model predictions using thermodynamic-based approaches have been of limited success in providing accurate predictions [32,33], probably due to missing controlling factors such as NADH oxidation and the dynamics of hydrogen partial pressure [31]. Another key factor explaining inaccurate predictions in existing rumen models is the lack of a dynamic representation of the microbial methanogens group. In this respect, new knowledge on the extent of methanogenesis and metabolic differences between the microbial members of this group could help to improve existing gut models. To our knowledge, none of the existing gut models integrates thermodynamics aspects to describe microbial growth of methanogens. Accordingly, our objective in this work was to develop a dynamic model with thermodynamic basis to quantify metabolism kinetics and energetic patterns that could inform on metabolic specificities between gut methanogenic archaea. Additionally, this model development aimed at providing tools for the analysis of ecological aspects (e.g., competitive exclusion principle) of the methanogenic community that can be
instrumental when designing, for example, nutritional strategies for methane mitigation in ruminants. The model was built upon quantitative data characterizing the *in vitro* dynamics of hydrogen utilization, methane production, growth and heat flux of three hydrogenotrophic methanogenic species representing major human and ruminant genera: *Methanobrevibacter smithii*, *Methanobrevibacter ruminantium* and *Methanobacterium formicium*.

# Material and methods

## *In vitro* growth experiments

**Archaeal strains and growth media.** Archaeal strains used in the study were *M. ruminantium* M1 (DSM 1093), *M. smithii* PS (type strain DSM 861), and *M. formicium* MF (type strain DSM 1535). All archaeal strains were purchased from DSM previously. The strains were frozen at -80°C, thawed and grown to get actively growing population before the beginning of the study. Balch growth media was prepared as previously described [34] and composition is summarized in S1 Table. Before autoclaving, 6 ml of media were distributed in Balch tubes (26 ml total volume) under CO₂ atmosphere.

**Experimental design and measures.** Starter cultures were grown until reaching optical density at 660 nm (OD₆₆₀) of 0.400 ± 0.030. Optical density was measured on a Jenway spectrophotometer (Bibby Scientific). Then, 0.6 ml were used to inoculate one experimental tube. Commercially prepared high purity H₂/CO₂ (80%/20%) gas mix was added to inoculated tubes by flushing for 1 min at 2.5 Pa. Mean initial OD₆₆₀ and pressure values are summarized in S2 Table. Growth kinetics for each strain were followed over 72 h. The experiment was repeated twice. Each kinetics study started with 40 tubes inoculated at the same time. At a given time point, two tubes with similar OD₆₆₀ values were sampled. The tubes were used for measuring gas parameters: pressure was measured using a manometer and composition of the gas phase was analysed by gas chromatography on a Micro GC 3000A (Agilent Technologies, France). The GC was equipped with two columns, MS-5A using argon as carrier gas and set to 100°C and PPU using helium and set to 75°C. The GC was calibrated using a certified gas standard mixture (Messer, France) containing methane, oxygen, hydrogen, carbon dioxide, and nitrogen. Approximately 2 ml of the sampled gas was injected in the GC for analysis. After gas sampling, one of the tubes was centrifuged 10 min at 13 000 g. The microbial pellet was weighed and stored at -20°C in 2 ml screw-cap tubes containing 0.4 g of sterile zirconia beads (0.3 g of 1 mm and 0.1 g of 0.5 mm).

**DNA extraction and qPCR quantification of 16S rRNA genes.** One ml of lysis buffer (50mM NaCl, 50 mM TrisHCl pH 7.6, 50 mM EDTA, 5% SDS) was added directly to the frozen microbial pellet before homogenizing for 2 × 30 s at 5100 tours/min in a Precellys bead-beater (Bertin Instruments). Samples were centrifuged for 3 min at 14 000 g and the liquid phase transferred to a new tube before adding 600 µl of phenol-chloroform-3-methyl-1-butanol (25:24:1) solution. After centrifugation at 14 000 g for 3 min, the aqueous phase was transferred to a fresh tube and 500 µl of chloroform were added. The chloroform-washing step was repeated twice with centrifugation at 14000 g for 3 min between steps. The final volume of the aqueous phase was measured and DNA precipitation was initiated by adding 70% of the volume of isopropanol 100% and 10% of the volume of sodium acetate 3M. Sedimentation at 14 000 g for 30 min was again performed and the resulting DNA pellet was washed with 500 µl of ethanol 70% and dissolved in 50µl of molecular biology quality water. The extraction yield was checked on a Nanodrop 1000 Spectrophotometer (Thermo Fisher Scientific, France) and extracts run on a FlashGel System (Lonza, Rockland, Inc) to check integrity.

Copies of 16S rRNA genes were quantified using a qPCR approach. Primers used are those of Ohene-Adjei et al [35]; reaction assay and temperature cycles were as described previously
[36]. Triplicate qPCR quantification was performed on 20 ng of extracted DNA. Amplifications were carried out using SYBR Premix Ex Taq (TaKaRa Bio Inc., Otsu, Japan) on a Ste- pOne system (Applied Biosystems, Courtabeuf, France). Absolute quantification involved the use of standard curves that had been prepared with gDNA of *Methanobrevibacter ruminant- ium* DSM 1093. PCR efficiency was of 103%. Results were expressed as copy numbers per ng of extracted DNA per g of microbial pellet. *M. smithii* and *M. ruminantium* strains used in this study possess two copies of 16S rRNA genes in their genomes. The number of cells was com- puted by dividing 16S copy numbers by 2.

Microcalorimetry

Microcalorimetric experiments were performed to determine the heat flux pattern of each methanogen. Metabolic activity and microbial growth were monitored by using isothermal calorimeters of the heat-conduction type (A TAM III, TA Instruments, France) equipped with two multicalorimeters, each holding six independent minicalorimeters, allowed contin- uous and simultaneous recording as a function of time of the heat flux produced by 12 sam- ples. The bath temperature was set at 39°C; its long-term stability was better than ± 1x10⁻⁴°C over 24h. Each minicalorimeter was electrically calibrated. The specific disposable 4 mL microcalorimetric glass ampoules capped with butyl rubber stoppers and sealed with alumin- ium crimps were filled with 1.75 mL of Balch growth media and overpressed with 2.5 Pa of H₂/CO₂ 80%/20% gas mixture for 30 s. There was no significant difference in pressure at the beginning of the study. They were sterilized by autoclave and stored at 39°C until the begin- ning of the microcalorimetric measurements. Actively growing cultures of methanogens (OD₆₆₀ of 0.280±0.030 for *M. smithii*, 0.271±0.078 for *M. ruminantium* and 0.142±0.042 for *M. formicium*) were stored at -20°C in order to diminish microbial activity before inocula- tion. Cultures were thawed for 30 min at ambient temperature and inoculation was carried out by injecting 0.25 mL of the culture through the septum of the overpressed microcalori- metric ampoules just before inserting them into the minicalorimeters. Samples took about two hours to reach the bath temperature and yield a stable zero baseline. Blank experiments were also carried out by inserting ampoules that were not inoculated and, as expected, no heat flux was observed confirming the medium sterility. Each experiment was repeated thrice.

The heat flux (dQ/dt), also called thermal power output P, was measured for each methanogen and blank samples with a precision ≥ 0.2 µW. The heat flux data of each sample were collected every 5 minutes during more than 10 days. The total heat Q was obtained by integrating the overall heat flux time curve using the TAM Assistant Software and its integrating function (TA Instruments, France).

Classically, the heat flux-time curve for a growing culture starts like the S-shaped biomass curve (a lag phase followed by an exponential growth phase) but differs beyond the growth phase, the heat flux being then modulated by transition periods [37]. Heat flux data can be used to infer the microbial growth rate constant provided the existence of a correlation between isothermal microcalorimetry data and microbiological data (e.g., cell counts) at early growth [38]. During the exponential growth phase, microbial growth follows a first-order kinetics defined by the specific growth rate constant µc (h⁻¹). Analogously, the heat flux follows an exponential behaviour determined by the parameter µc as described by [37,38].

$$ \frac{dQ}{dt} = \mu_c \cdot Q \qquad (1) $$

The growth rate constant µc can be determined by fitting the exponential part of the
<figure><img src="image_5.png" /><figcaption>**Fig 1.** Schematics of the *in vitro* methanogenesis process by hydrogenotrophic methanogens. Double arrows represent fluxes due to liquid-gas transfer, simple arrows represent metabolic fluxes.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0226243.g001

heat flux-time curve using the fitting function of the TAM Assistant Software. In our case study, careful selection of the exponential phase of heat flux dynamics was performed to provide a reliable estimation of the maximum growth rate constant from calorimetric data.

**Mathematical model development**

**Modelling in vitro methanogenesis.** The process of *in vitro* methanogenesis is depicted in Fig 1. The H₂/CO₂ mixture in the gas phase diffuses to the liquid phase. The H₂ and CO₂ in the liquid phase are further utilized by the pure culture to produce CH₄. Methane in the liquid phase diffuses to the gas phase.

Model construction was inspired by our previous dynamic models of human gut [19] and rumen *in vitro* fermentation [39] followed by certain simplifications. The model considers the liquid-gas transfer of carbon dioxide. Due to the low solubility of hydrogen and methane [40], the concentration of these two gases in the liquid phase was not modelled. We assumed that the dynamics of concentrations in the gas phase are determined by kinetic rate of the methanogenesis. To incorporate thermodynamic information, instead of using the Monod equation in the original formulation, we used the kinetic rate function proposed by Desmond-Le Quéméner and Bouchez [26]. The resulting model is described by the following ordinary differential
equations

$$ \frac{dx_{H_2}}{dt} = \mu_{\max} \cdot \exp \left( - \frac{K_s \cdot V_g}{n_{g,H_2}} \right) \cdot x_{H_2} - k_d \cdot x_{H_2} \quad (2) $$

$$ \frac{ds_{CO_2}}{dt} = - \frac{Y_{CO_2} \mu_{\max}}{Y} \cdot \exp \left( - \frac{K_s \cdot V_g}{n_{g,H_2}} \right) \cdot x_{H_2} - k_L a \cdot \left( s_{CO_2} - K_{H,CO_2} \cdot R \cdot T \cdot n_{g,CO_2} / V_g \right) \quad (3) $$

$$ \frac{dn_{g,H_2}}{dt} = - \frac{\mu_{\max}}{Y} \cdot \exp \left( - \frac{K_s \cdot V_g}{n_{g,H_2}} \right) \cdot V_L \cdot x_{H_2} \quad (4) $$

$$ \frac{dn_{g,CO_2}}{dt} = V_L \cdot k_L a \cdot \left( s_{CO_2} - K_{H,CO_2} \cdot R \cdot T \cdot n_{g,CO_2} / V_g \right) \quad (5) $$

$$ \frac{dn_{g,CH_4}}{dt} = \frac{Y_{CH_4} \mu_{\max}}{Y} \cdot \exp \left( - \frac{K_s \cdot V_g}{n_{g,H_2}} \right) \cdot V_L \cdot x_{H_2} \quad (6) $$

where $s_{CO_2}$ is the concentration (mol/L) of carbon dioxide in the liquid phase and $x_{H_2}$ is the biomass concentration (mol/L) of hydrogenotrophic methanogens. The numbers of moles in the gas phase are represented by the variables $n_{g,H_2}$, $n_{g,CO_2}$, $n_{g,CH_4}$. The gas phase volume $V_g = 20$ mL and the liquid phase volume $V_L = 6$ mL. Liquid-gas transfer for carbon dioxide is described by a non-equilibria transfer rate which is driven by the gradient of the concentration of the gases in the liquid and gas phase. The transfer rate is determined by the mass transfer coefficient $k_L a (h^{-1})$ and the Henry's law coefficients $K_{H,CO_2}$ (M/bar). $R (bar \cdot (M \cdot K)^{-1})$ is the ideal gas law constant and $T$ is the temperature (K). Microbial decay is represented by a first-order kinetic rate with $k_d (h^{-1})$ the death cell rate constant. Microbial growth was represented by the rate function proposed by Desmond-Le Quéméner and Bouchez [26] using hydrogen as single substrate

$$ \mu = \mu_{\max} \cdot \exp \left( - \frac{K_s \cdot V_g}{n_{g,H_2}} \right) \quad (7) $$

where $\mu$ is the growth rate (h$^{-1}$), $\mu_{\max}$ (h$^{-1}$) is the maximum specific growth rate constant and $K_s (mol/L)$ the affinity constant. Eq (7) is derived from energetic principles following Boltzmann statistics and uses the concept of exergy (maximum work available for a microorganism during a chemical transformation). The affinity constant has an energetic interpretation since it is defined as

$$ K_s = \frac{E_{M+} E_{dis}}{v_{harv} \cdot E_{cat}} \quad (8) $$

where $E_{dis}$ (kJ/mol) and $E_M$ (kJ/mol) are the dissipated exergy and stored exergy during growth respectively. $E_{cat}$ (kJ/mol) is the catabolic exergy of one molecule of energy-limiting substrate, and $v_{harv}$ is the volume at which the microbe can harvest the chemical energy in the form of substrate molecules [26]. $E_{cat}$ is the absolute value of the Gibbs energy of catabolism ($\Delta G_{r,c}$) when the reaction is exergonic ($\Delta G_{r,c} < 0$) or zero otherwise. The stored exergy $E_M$ is calculated from a reaction (destock) representing the situation where the microbe gets the energy by consuming its own biomass. $E_M$ is the absolute value of the Gibbs energy of biomass consuming
reaction (ΔG<sub>r,destock</sub>) when the reaction is exergonic (ΔG<sub>r,destock</sub> < 0) or zero otherwise. Finally, the dissipated exergy E<sub>dis</sub> is the opposite of the Gibbs energy of the overall metabolic reaction, which is a linear combination of the catabolic and destock reactions. This calculation follows the Gibbs energy dissipation detailed in Kleerebezem and Van Loosdrecht [24].

In our model, the stoichiometry of methanogenesis is represented macroscopically by one catabolic reaction (R<sub>1</sub>) for methane production and one anabolic reaction (R<sub>2</sub>) for microbial formation. We assumed that ammonia is the only nitrogen source for microbial formation. The molecular formula of microbial biomass was assumed to be C<sub>5</sub>H<sub>7</sub>O<sub>2</sub>N [40].

$$R_1 : 4 H_2 + CO_2 \rightarrow CH_4 + 2 H_2O$$

$$R_2 : 10 H_2 + 5 CO_2 + NH_3 \rightarrow C_5H_7O_2N + 8 H_2O$$

In the model, the stoichiometry of the reactions is taken into account via the parameters Y, Y<sub>CO₂</sub>, Y<sub>CH₄</sub>, which are the yield factors (mol/mol) of microbial biomass, CO₂ and CH₄. The microbial yield factor Y was extracted from literature. We assumed that M. smithii and M. ruminantium have the same yield (being both Methanobrevibacter). This yield factor was set to 0.006 mol biomass/mol H₂, using the methane-based molar growth yield of 2.8 g biomass/mol CH₄ estimated for M. smithii [41] and the Eqs (9) and (11). Similarly, the yield factor for M. formicium was set to 0.007 mol biomass/mol H₂ using the methane-based molar growth yield of 3.5 g biomass/mol CH₄ reported by Schauer and Ferry [42]. The fraction of H₂ utilized for microbial growth (reaction R₂) is defined by the yield factor Y (mol of microbial biomass/mol of H₂). Now, let f be the fraction of H₂ used for the catabolic reaction R₁. Reaction R₂ tells us that for every 10 moles of H₂ used in R₂, we get 1 mol of microbial biomass. Hence, it follows that

$$Y = \frac{1}{10} \cdot (1 - f) \qquad (9)$$

If Y is known, the fraction f can be computed from Eq (9).

The yield factors of CO₂ and CH₄ can be expressed as functions of the fraction f:

$$Y_{CO_2} = \left(\frac{1}{4}\right) \cdot f + \left(\frac{5}{10}\right) \cdot (1 - f) \qquad (10)$$

$$Y_{CH_4} = \left(\frac{1}{4}\right) \cdot f \qquad (11)$$

The model has two physicochemical parameters (k<sub>L</sub>a, K<sub>H,CO₂</sub>) and four biological parameters (μ<sub>max</sub>, K<sub>s</sub>, Y, k<sub>d</sub>). The initial condition for s<sub>CO₂</sub> is unknown and was also included in the parameter vector for estimation. The Henry's law coefficients are known values calculated at 39°C using the equations provided by Batstone et al. [40]. An implementation of the model in the Open Source software Scilab (https://www.scilab.org/) is available at https://doi.org/10.5281/zenodo.3271611.

**Theoretical model to study interactions among methanogens.** The experimental study of microbial interactions requires sophisticated *in vitro* systems under continuous operation such as the one developed by Haydock *et al.* [43]. In our work, we explored by means of mathematical modelling how the methanogens can interact under *in vivo* conditions. For this theoretical study, we elaborated a toy model based on the previous model for *in vitro* methanogenesis. Let us consider the following simple model for representing the consumption of
hydrogen by the methanogenic species i under an in vivo scenario of continuous flow

$$ \frac{dx_{H_{2,i}}}{dt} = \mu_{\max,i} \cdot \exp\left(-\frac{K_{s,i} \cdot V_g}{n_{g,H_2}}\right) \cdot x_{H_{2,i}} - D_i \cdot x_{H_{2,i}} \quad (12) $$

$$ \frac{dn_{g,H_2}}{dt} = q_{H_2} - \frac{\mu_{\max,i}}{Y_i} \cdot \exp\left(-\frac{K_{s,i} \cdot V_g}{n_{g,H_2}}\right) \cdot V_L \cdot x_{H_2,i} - b \cdot n_{g,H_2} \quad (13) $$

where $q_{H_2}$ (mol/h) is the flux of hydrogen produced from the fermentation of carbohydrates. The kinetic parameters are specific to the species i ($x_{H_2,i}$). The parameter $D_i$ (h⁻¹) is the dilution rate of the methanogens and b (h⁻¹) is an output substrate rate constant. Extending the model to n species with a common yield factor Y, the dynamics of hydrogen is given by

$$ \frac{dn_{g,H_2}}{dt} = q_{H_2} - \frac{V_L}{Y} \sum_{i=1}^{n} \mu_{\max,i} \cdot \exp\left(-\frac{K_{s,i} \cdot V_g}{n_{g,H_2}}\right) \cdot x_{H_2,i} - b \cdot n_{g,H_2} \quad (14) $$

where the sub index i indicates the species. In our case study, n = 3.

## Parameter identification

Before performing the numerical estimation of the model parameters, we addressed the question of whether it was theoretically possible to determine uniquely the model parameters given the available measurements from the experimental setup. This question is referred to as structural identifiability [44]. Structural identifiability analysis is of particular relevance for model whose parameters are biologically meaningful, since knowing the actual value of the parameter is useful for providing biological insight on the system under study [45]. Moreover, in our case, we are interested in finding accurate estimates that can be further used as priors in an extended model describing the *in vivo* system.

We used the freely available software DAISY [46] to assess the structural identifiability of our model. Physical parameters ($k_{L}$a, $K_{H,CO_2}$) were set to be known. The model was found to be structurally globally identifiable. In practice, however, to facilitate the actual identification of parameters and reduce practical identifiability problems such as high correlation between the parameters [47], we fixed some model parameters to values reported in the literature. The transport coefficient $k_{L}$a, the Henry's law coefficient $K_{H,CO_2}$, and the dead cell rate constant $k_{d}$ were set to be known and were extracted from Batstone et al. [40]. Therefore, only the parameters $\mu_{\max}$, $K_{s}$ and initial condition of $s_{CO_2}$ were set to be estimated. To capitalize on the calorimetric data, we further assumed that $\mu_{\max}$ was equal to the specific rate constant $\mu_{c}$ estimated from the heat flux-time curve. By this, only the affinity constant for each strain and the initial condition of $s_{CO_2}$ were left to be estimated.

The parameter identification for each methanogen was performed with the IDEAS Matlab® toolbox [48] (freely available at http://genome.jouy.inra.fr/logiciels/IDEAS). The parameter identification was performed with the data of the *in vitro* growth experiments. The measured variables are the number of moles in the gas phase (H₂, CH₄, CO₂). The Lin's concordance correlation coefficient (CCC) [49] was computed to quantify the agreement between the observations and model predictions.
# Results

## Methanogens biomass

Archaea-specific primers targeting the 16S rRNA gene were used to enumerate microbial cells in each pure culture. Three hours post inoculation microbial numbers varied from $7.62×10^7$ to $2.81×10^8$ and reached $10^9$ after 72 hours of incubation. S3 Table summarizes microbial numbers at different sampling times.

## Calorimetric pattern of methanogens

Fig 2 displays a representative isothermal calorimetric curve for each methanogen. The three measured heat flux dynamics of each methanogen were found to follow similar energetic patterns. *M. smithii* and *M. formicium* exhibited a lag phase of a few hours, while *M. ruminantium* was already metabolically active when introduced into the minicalorimeter though several attempts were made to obtain a lag phase by changing storage conditions and thawing the

<figure><img src="image_6.png" /><figcaption>**Fig 2.** Example of isothermal calorimetric curves for *M. ruminantium* (dashed blue line), *M. smithii* (solid red line) and *M. formicium* (dotted black line). The dominant metabolic phase is represented by one peak (shown with the arrows). The magnitude of the peak differs between the methanogens and also the slope of the heat flux trajectories. The return of the heat flux to the zero baseline also differs between the three methanogens. The inset zoom displays the peak exhibited by *M. smithii* at 7.4 h.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0226243.g002
culture just before inoculating the microcalorimetric ampoules. The pattern of heat flux for all tested methanogens is characterized by one predominant peak which was observed at different times for each methanogen. *M. smithii* exhibited a second metabolic event occurring at 60 h with an increase of heat flux. The same phenomenon was observed for *M. formicium* but at a lower intensity that started at 140 h. The process was considered completed when the heat flux ceased marking the end of the metabolic activity. It is noted that *M. formicium* produced a small peak at 14 h (Fig 2). A similar peak, but of much smaller size, was observed on the other curves obtained with this methanogen. *M. smithii* also exhibits a small peak (occurrence of 3 out of 3) at 7.4 h shown in the inset of Fig 2.

The total heat ($Q_m$) produced during the methanogenesis process that took place under the present experimental conditions was, on average, -5.5 ± 0.5 J for the three methanogens (for *M. ruminantium*, the missing initial part of the heat flux-time curve was approximately estimated by extrapolating the exponential fit). As we shall see below, this experimental value is consistent with the theoretically expected value.

## Estimation of thermodynamic properties

We defined two macroscopic reactions to represent the catabolism (R1) and anabolism (R2) of the methanogenesis (see Modelling in vitro methanogenesis section). All thermodynamic properties result from the contribution of both catabolic and anabolic reactions. The calculations of total heat ($Q_m$), enthalpy ($\Delta H_m$), Gibbs energy ($\Delta G_m$) and entropy ($\Delta S_m$) of the methanogenesis are detailed in S4 Table. The estimated overall heat produced during the methanogenesis process under our experimental conditions was in average $Q_m = -5.66$ J. This heat results from the sum of the heat of the catabolic reaction ($Q_c$) and the heat of the anabolic reaction ($Q_a$). From the total heat of the methanogenesis, the anabolic reaction contributes to 7% of the metabolic heat for *M. smithii* and *M. ruminantium*. For *M. formicium*, the contribution of the anabolic reaction to the metabolic heat is 9%. It is also interesting to note that there is a very good agreement between the theoretical value calculated above and the overall heat experimentally determined by microcalorimetry (-5.5 ± 0.5 J).

Table 1 shows the thermodynamic properties per mole of biomass formed during methanogenesis of *M. ruminantium*, *M. smithii* and *M. formicium* on $H_2/CO_2$. These properties are compared with values found in the literature for other methanogens grown on different substrates.

## Dynamic description of *in vitro* kinetics

The developed mathematical model was calibrated with the experimental data from *in vitro* growth experiments in Balch tubes. Table 2 shows the parameters of the dynamic kinetic model described in Eqs 2–6. The reported value of $\mu_{max}$ for each methanogen corresponds to

<table><caption>Table 1. Gibbs energies, enthalpies and entropies of metabolic processes involving some methanogens growing on different energy sources.</caption><thead><tr><th>Microorganism</th><th>Energy substrate</th><th>Growth conditions</th><th>ΔG<sub>m</sub><br>kJ / C-mol</th><th>ΔH<sub>m</sub><br>kJ / C-mol</th><th>TΔS<sub>m</sub><br>kJ / K<sup>-1</sup> C-mol</th><th>Driving force</th><th>Reference</th></tr></thead><tbody><tr><td><i>M. ruminantium</i>, <i>M. smithii</i></td><td>H<sub>2</sub>/CO<sub>2</sub></td><td>anaerobic</td><td>-1073</td><td>-2132</td><td>-1059</td><td>Enthalpy-driven but Entropy-retarded</td><td>this work</td></tr><tr><td><i>M. formicium</i></td><td>H<sub>2</sub>/CO<sub>2</sub></td><td>anaerobic</td><td>-801</td><td>-1605</td><td>-804</td><td>Enthalpy-driven but Entropy-retarded</td><td>this work</td></tr><tr><td><i>M. thermo-autotrophicum</i></td><td>H<sub>2</sub>/CO<sub>2</sub></td><td>anaerobic</td><td>-802</td><td>-3730</td><td>-2928</td><td>Enthalpy-driven but Entropy-retarded</td><td>[50]</td></tr><tr><td><i>M. formicium</i></td><td>formate</td><td>anaerobic</td><td>-880</td><td>-613</td><td>+267</td><td>Enthalpy-driven</td><td>[23]</td></tr><tr><td><i>M. barkeri</i></td><td>methanol</td><td>anaerobic</td><td>-570</td><td>-420</td><td>+150</td><td>Enthalpy-driven</td><td>[23]</td></tr><tr><td><i>M. barkeri</i></td><td>acetate</td><td>anaerobic</td><td>-366</td><td>+145</td><td>+511</td><td>Entropy-driven but enthalpy-retarded</td><td>[51]</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0226243.t001
<table><caption>**Table 2. Parameters of the model of in vitro methanogenesis.** The value reported $\mu_{max}$ for each methanogen is the mean value obtained from three heat flux-time curves.</caption><thead><tr><th>Parameter</th><th>Definition</th><th colspan="3">Value</th></tr></thead><tbody><tr><td>$k_{1a}$ (h<sup>-1</sup>)</td><td>Liquid-gas transfer constant</td><td colspan="3">8.33</td></tr><tr><td>$K_{H,CO_2}$ (M/bar)</td><td>Henry's law coefficient of carbon dioxide</td><td colspan="3">0.0246</td></tr><tr><td>$k_d$ (h<sup>-1</sup>)</td><td>Death cell rate constant</td><td colspan="3">8.33x10<sup>-4</sup></td></tr><tr><td></td><td></td><th>M. smithii</th><th>M. ruminantium</th><th>M. formicium</th></tr><tr><td>$K_s$ (mol/L)</td><td>Affinity constant</td><td>0.028</td><td>0.042</td><td>0.011</td></tr><tr><td>$\mu_{max}$ (h<sup>-1</sup>)</td><td>Maximum specific growth rate constant</td><td>0.12</td><td>0.07</td><td>0.046</td></tr><tr><td>Y (mol biomass/mol H<sub>2</sub>)</td><td>Microbial biomass yield factor</td><td>0.006</td><td>0.006</td><td>0.007</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0226243.t002

the average value obtained from three heat flux-time curves. From Table 2, it is concluded that *M. smithii* exhibited the highest growth rate constant, followed by *M. ruminantium* and finally *M. formicium*. In terms of the affinity constant $K_s$, while *M. smithii* and *M. ruminantium* are of the same order, the affinity constant for *M. formicium* is lower in one order of magnitude.

Fig 3 displays the dynamics of the compounds in the methanogenesis for the three methanogens. The data of the figure is available at https://doi.org/10.5281/zenodo.3469655. Experimental data are compared against the model responses. Table 3 shows standard statistics for model evaluation. The model captures efficiently the overall dynamics of the methanogenesis. Hydrogen and methane are very well described by the model with average concordance correlation coefficients (CCC) of 0.98 and 0.97 respectively. For carbon dioxide, CCC = 0.95.

Fig 4 displays the dynamics for the methanogens as measured by the 16S rRNA gene, as well as the dynamics of biomass as predicted for the model. As observed, the microbes follow a typical Monod-like trajectory.

<figure><img src="image_7.png" /><figcaption>**Fig 3. Assessment of model performance.** Top plots: dynamics of methanogenesis by *M. ruminantium* (*), *M. smithii* (○) and *M. formicium* (□). Experimental data (*,○,□) are compared against model predicted responses: dotted blue lines (*M. ruminantium*), solid red lines (*M. smithii*) and dashed black lines (*M. formicium*). Bottom plots: summary observed vs predicted variables. The solid red line is the isocline.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0226243.g003
<table><caption>Table 3. Statistical indicators for model evaluation.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Hydrogen</th><th colspan="3">Methane</th><th colspan="3">Carbon dioxide</th></tr><tr><th>CCC*</th><th>r<sup>2</sup></th><th>CV<sub>RMSE</sub><sup>**</sup></th><th>CCC*</th><th>r<sup>2</sup></th><th>CV<sub>RMSE</sub><sup>**</sup></th><th>CCC*</th><th>r<sup>2</sup></th><th>CV<sub>RMSE</sub><sup>**</sup></th></tr></thead><tbody><tr><td><i>M. smithii</i></td><td>0.99</td><td>0.98</td><td>16</td><td>0.96</td><td>0.94</td><td>18</td><td>0.93</td><td>0.84</td><td>6</td></tr><tr><td><i>M. ruminantium</i></td><td>0.97</td><td>0.95</td><td>11</td><td>0.96</td><td>0.93</td><td>20</td><td>0.99</td><td>0.98</td><td>2</td></tr><tr><td><i>M. formicium</i></td><td>0.99</td><td>0.98</td><td>13</td><td>0.98</td><td>0.96</td><td>17</td><td>0.92</td><td>0.86</td><td>8</td></tr><tr><td>Mean</td><td>0.98</td><td>0.97</td><td>13</td><td>0.97</td><td>0.94</td><td>18</td><td>0.95</td><td>0.89</td><td>5</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0226243.t003

## Discussion

Our objective in this work was to quantitatively characterize the dynamics of hydrogen utilization, methane production, growth and heat flux of three hydrogenotrophic methanogens by integrating microbiology, thermodynamics and mathematical modelling. Our model developments were instrumental to quantify energetic and kinetic differences between the three methanogens studied, strengthening the potentiality of microcalorimetry as a tool for characterizing the metabolism of microorganisms [52]. This modelling work provides estimated parameters that can be used as prior values for other modelling developments of gut microbiota.

## Energetic and kinetic differences between methanogens

Methanogenesis appears as a simple reaction described by a single substrate kinetic rate on H₂. The microcalorimetry approach we applied revealed that this simplicity is only apparent and that hydrogenotrophic methanogens exhibit energetic and kinetic differences. Methanogenesis is indeed a complex process that can be broken down in several stages. The dominant metabolic phase is represented by one peak that occurs at different times. The magnitude of the peak differs between the methanogens and also the slope of the heat flux trajectories. For *M. smithii* and *M. formicium* the main peak was preceded by a small increase in heat flux which translates in a metabolic activity that remains to be elucidated. For *M. ruminantium*, we do not

<figure><img src="image_8.png" /><figcaption>Fig 4. Dynamics of methanogens as measured by 16S rRNA gene copies (circles) and biomass concentrations (solid line) predicted by the model.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0226243.g004

<footnote>* CCC: Lin's concordance correlation coefficient.</footnote>

<footnote>** CV<sub>RMSE</sub>: coefficient of variation of the root mean squared error.</footnote>
know whether the small peak exists since the initial part of the curve is missing. For *M. smithii* and *M. formicium*, it was observed a metabolic event after the main peak around 60 h for *M. smithii* and 140h *M. formicium*. This peak could be representing cell lysis process [38]. The return time of the heat flux to the zero baseline was also different. The energetic difference is associated with kinetic differences that translate into specific kinetic parameters, namely affinity constant ($K_s$) and maximum growth rate constant ($\mu_{max}$). Previously, energetic differences between methanogens have been ascribed to the presence or absence of cytochromes [15]. These differences are translated into different yield factors, H₂ thresholds, and doubling times. The kinetic differences revealed in this study for three cytochrome-lacking methanogens indicate that factors other than the presence of cytochromes might play a role in the energetics of methanogenesis. Interestingly, calorimetric experiments showed that *M. ruminantium* was metabolically active faster than the other methanogens, characteristic that could explain the predominance of *M. ruminantium* in the rumen [53]. Looking at the expression of the affinity constant (Eq (8)), the differences between the affinity constants among the methanogens can be explained by the differences between the by the harvest volume $v_{harv}$ and the yield factors. Note that in the kinetic function developed by Desmond-Le Quéméner and Bouchez [26], the maximum growth rate did not have any dependency on the energetics of the reaction. Our experimental study revealed that $\mu_{max}$ is species-specific and reflects the dynamics of the heat flux of the reaction at the exponential phase. This finding suggests that a further extension of the kinetic model developed by Desmond-Le Quéméner and Bouchez [26] should include the impact of energetics on $\mu_{max}$. Since our study is limited to three species, it is important to conduct further research on other methanogens to validate our findings. In this same line, to enhance the evaluation of the predictive capabilities of our model, a further model validation is required with independent data set under different experimental conditions (e.g. continuous mode operation) to those used in this study.

## Energetic analysis

Regarding the energetic information for different methanogens summarized in Table 1, it is observed that the thermodynamic behaviour of the three methanogens is analogous to that observed for *Methanobacterium thermoautotrophicum* [50]. The values reported in Table 1 show indeed that the methanogenesis on H₂/CO₂ is characterized by large heat production. The growth is highly exothermic, with a $\Delta H_m$ value that largely exceeds the values found when other energy substrates are used. The enthalpy change $\Delta H_m$, which is more negative than the Gibbs energy change $\Delta G_m$, largely controls the process. Growth on H₂/CO₂ is also characterized by a negative entropic contribution $T\Delta S_m$ which, at first sight, may look surprising since entropy increases in most cases of anaerobic growth [54]. However, this can be understood if one remembers that $T\Delta S_m$ corresponds in fact to the balance between the final state and the initial state of the process, that is

$$T\Delta S_m = \frac{(1 - 10Y)}{4Y} T\Delta S_c + T\Delta S_a = \frac{(1 - 10Y)}{4Y} T(S_{final} - S_{initial})_c + T(S_{final} - S_{initial})_a$$

Methanogenesis on H₂/CO₂ is particular because the final state of its catabolic reaction (1 mol CH₄ + 2 mol H₂O) involves a smaller number of moles than the initial state (4 mol H₂ + 1 mol CO₂), which results in a significant loss of entropy during the process. For spontaneous growth in such a case, the $\Delta H_m$ must not only contribute to the driving force but must also compensate the growth-unfavourable $T\Delta S_m$, which means that $\Delta H_m$ must be much more negative than $\Delta G_m$ [55]. For this reason, methanogenesis on H₂/CO₂, which is accompanied by a considerable decrease of entropy and a large production of heat, has been designed as an
entropy-retarded process [50]. More generally, von Stockar and Liu [55] noticed that when the Gibbs energy of the metabolic process is resolved into its enthalpic and entropic contributions, very different thermodynamic behaviours are observed depending on the growth type. These thermodynamic behaviours are: aerobic respiration is clearly enthalpy-driven (ΔHm ≪ 0 and TΔSm > 0), whereas fermentative metabolism is mainly entropy-driven (ΔHm < 0 and TΔSm ≫ 0). Methanogenesis on H₂/CO₂ is enthalpy-driven but entropy-retarded (ΔHm ≪ 0 and TΔSm < 0), whereas methanogenesis on acetate is entropy-driven but enthalpy-retarded (ΔHm > 0 and TΔSm ≫ 0). In the present case, the highly exothermic growth of M. ruminantium, M. smithii and M. formicium on H₂/CO₂ is largely due to the considerable decrease of entropy during the process: in fact, 50% of the heat produced here serves only to compensate the loss of entropy. A proportion of 80% was found for M. thermoautotrophicum [50], which results from the fact that their TΔSm and ΔHm values are, respectively, 2.7 and 1.7 times larger than ours. This difference might be due to the differences in the temperature of the studies, namely 39°C in our study vs 60°C in the study by Schill et al. [50].

## Do our results inform on ecological questions such as species coexistence?

The competitive exclusion principle [56] states that coexistence cannot occur between species that occupy the same niche (the same function). Only the most competitive species will survive. Recently, by using thermodynamic principles, Großkopf & Soyer [27] demonstrated theoretically that species utilizing the same substrate and producing different compounds can coexist by the action of thermodynamic driving forces. Since in our study the three methanogens perform the same metabolic reactions, the thermodynamic framework developed Großkopf & Soyer [27] predicts, as the original exclusion principle [56], the survival of only one species. By incorporating thermodynamic control on microbial growth kinetics, Lynch et al [57] showed theoretically that differentiation of ATP yields can explain ecological differentiation of methanogens over a range of liquid turnover rates. This theoretical work predicts that for a fixed liquid turnover rate, only one species survives. For the continuous culture of microorganisms, it has been demonstrated that at the equilibrium (growth rate equals the dilution rate) with constant dilution rates and substrate input rates, the species that has the lowest limiting substrate concentration wins the competition. From Eq (12), the number of moles of hydrogen of the species $n_{g,H_2,i}^*$ at the steady state is

$$n_{g,H_2,i}^* = \frac{K_{s,i} \cdot V_g}{\log(\mu_{max,i}/D_i)}$$

Using the model parameters of Table 2, we studied in silico three possible competition scenarios, assuming a constant environment (constant dilution rate D). Two dilution rates were evaluated: D = 0.021 h⁻¹ (retention time = 48 h) and D = 0.04 h⁻¹ (retention time = 25 h). A retention time of 48 h corresponds to values measured in small ruminants [58] and to humans as we used in our gut model [19]. For higher retention times, the results obtained for 48 h hold. For D = 0.021 h⁻¹, we obtained that $n_{g,H_2,Ms}^* = 0.32$ mmol, $n_{g,H_2,Mr}^* = 0.68$ mmol, $n_{g,H_2,Mf}^* = 0.28$ mmol, where the subindex Ms, Mr, Mf stand for M. smithii, M. ruminantium and M. formicium. From these results, it appears that under a constant environment, M. formicium will win the competition. Since $n_{g,H_2,Ms}^* < n_{g,H_2,Mr}^*$, M. ruminantium will be extinguished before M. smithii. For D = 0.04 h⁻¹, we obtained that $n_{g,H_2,Ms}^* = 0.49$ mmol, $n_{g,H_2,Mr}^* = 1.42$ mmol, $n_{g,H_2,Mf}^* = 1.57$ mmol, and thus M. smithii wins the competition. To win the competition, M. ruminantium requires longer retention times than its competitors. Retention times of digesta longer than 48 h are physiologically uncommon, thus the presence
<figure><img src="image_9.png" /><figcaption>**Fig 5.** Possible competition scenarios between *M. ruminantium* (blue dashed line), *M. smithii* (red solid line) and *M. formicium* (black dotted line) in a hypothetical constant environment. A. At constant dilution rate of 0.021 h⁻¹, *M. formicium* displaces the other two methanogens. B. With a constant dilution rate of 0.04 h⁻¹, *M. smithii* wins the competition. At constant environmental conditions, only one species wins and displaces the other methanogens.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0226243.g005

of *M. ruminantium* in the gut ecosystem can be explained, for example, from known adhesion properties (both *M. ruminantium* and *M. smithii* genes encode adhesin-like proteins [59,60]). To illustrate these aspects, we built a multiple-species model with the three methanogens using Eqs (12) and (14). The parameter *b* was set to 0.5 h⁻¹ and the hydrogen flux production *q*<sub>H₂</sub> rate was set to 0.02 mol/min. Fig 5A displays the dynamics of the three methanogens for the first scenario (D = 0.021 h⁻¹). It is observed that at 50 d only *M. formicium* survives. This result, however, is not representative of what occurs in the rumen where the three methanogens coexist [5,61]. It is intriguing that in our toy model it is *M. formicium* that wins the competition, bearing in mind that *M. ruminantium* and *M. smithii* are more abundant than *M. formicium* [5,53]. Fig 5 shows that selective conditions favour the survival of one species. Similar results can be obtained for the human gut by including the effect of pH on microbial growth [22] and setting the gut pH to select one of the species.

On the basis of the competitive exclusion principle, it is thus intriguing that having a very
specialized function, methanogens are a diverse group that coexist. Gut ecosystems, therefore,
exhibit the paradox of the plankton introduced by Hutchinson (1961) that presents the coexis-
tence of species all competing for the same substrate in a relatively isotropic or unstructured
environment [62]. In the case of the rumen, our modelling work suggests that in addition to
kinetic and thermodynamic factors, other forces contribute to the ecological shaping of the
methanogens community in the rumen favouring the microbial diversity. Indeed, methano-
genic diversity in the rumen results from multiple factors that include pH sensitivity, the asso-
ciation with rumen fractions (fluid and particulate material), and the endosymbiosis with
rumen protozoa [5,53]. For the human gut, ecological factors enable methanogens to coexist
to a competitive environment where hydrogenotrophic microbes (acetogens, methanogenic
archaea and sulfate-reducing bacteria) utilize H₂ *via* different pathways [63–65]. Both in the human gut and in the rumen, microbes grow in association with biofilms that form a polymer-based matrix that provides nutritional and hydraulic advantages for microbial growth and resistance to shear forces [19,66]. Indeed, in our modelling work of human gut fermentation [19], we suggested that, from the different actions the mucus has on colonic fermentation, the mechanism of promoting conditions for microbial aggregation appears as the most relevant factor for attaining the high microbial density and the high level of fibre degradation characteristic of the human gut. Altogether, these factors result in nonlinear behaviours, spatial and temporal variations that promote coexistence and diversity, that, as discussed in dedicated literature on microbial ecology [67–71], render the classical formulation of the competitive exclusion principle [56,72] inapplicable to gut ecosystems.

Finally, mathematical modelling is expected to enhance our understanding of gut ecosystems [66,73]. It is then key that in addition to metabolic aspects, mathematical models of gut fermentation incorporate the multiple aspects that shape microbial dynamics to provide accurate predictions and improve insight on gut metabolism dynamics and its potential modulation. For ruminants, the development of precision livestock technologies provides promising alternatives for integrating real-time data of key animal phenotypes such as feeding behaviour with mathematical models for estimating methane emissions [74] and rumen function indicators at large scale. These tools will be instrumental to support livestock management decisions and guide timely interventions. Similarly, for humans, mathematical models coupled with electronic technologies for online monitoring of gut function [75] might facilitate the diagnosis and the design of personalized therapies for gastrointestinal diseases.

**Supporting information**

**S1 Table. Methanogens growth media composition.**

(DOCX)

**S2 Table. Summary of initial OD and pressure measured immediately after primary inoculation.**

(DOCX)

**S3 Table. qPCR quantification of 16S rRNA genes.**

(DOCX)

**S4 Table. Calculation of thermodynamic properties of the methanogenesis.**

(DOCX)

**Acknowledgments**

We are grateful to Dominique Graviou (UMRH, Inra) for her skilled assistance on the *in vitro* growth experiments and qPCR assays. We thank the Inra PHASE department and the Inra MEM metaprogramme for financial support. RMT, MP and DPM acknowledge the support of ERA-net gas co-fund for funding the RumenPredict project.

**Author Contributions**

**Conceptualization:** Rafael Muñoz-Tamayo, Milka Popova, Diego P. Morgavi, Nicole Morel-Desrosiers.

**Data curation:** Milka Popova, Nicole Morel-Desrosiers.

**Formal analysis:** Rafael Muñoz-Tamayo, Milka Popova, Nicole Morel-Desrosiers.
**Funding acquisition:** Rafael Muñoz-Tamayo, Milka Popova.

**Investigation:** Rafael Muñoz-Tamayo, Milka Popova, Maxence Tillier, Nicole Morel-Desrosiers.

**Methodology:** Rafael Muñoz-Tamayo, Milka Popova, Maxence Tillier, Diego P. Morgavi, Jean-Pierre Morel, Gérard Fonty, Nicole Morel-Desrosiers.

**Software:** Rafael Muñoz-Tamayo.

**Writing – original draft:** Rafael Muñoz-Tamayo, Milka Popova, Nicole Morel-Desrosiers.

**Writing – review & editing:** Rafael Muñoz-Tamayo, Milka Popova, Diego P. Morgavi, Jean-Pierre Morel, Gérard Fonty, Nicole Morel-Desrosiers.

## References

1. Miller TL, Wolin MJ, Hongxue Z, Bryant MP. Characteristics of methanogens isolated from bovine rumen. Appl Environ Microbiol. 1986; 51: 201–202. PMID: 3954338

2. Dridi B, Fardeau ML, Ollivier B, Raoult D, Drancourt M. Methanomassiliicoccus luminyensis gen. nov., sp. nov., a methanogenic archaeon isolated from human faeces. Int J Syst Evol Microbiol. 2012; 62: 1902–1907. https://doi.org/10.1099/ijms.0.033712-0 PMID: 22859731

3. Paul K, Nonoh JO, Mikulski L, Brune A. "Methanoplasmatales," thermoplasmatales-related archaea in termite guts and other environments, are the seventh order of methanogens. Appl Environ Microbiol. 2012; 78: 8245–8253. https://doi.org/10.1128/AEM.02193-12 PMID: 23001661

4. Dridi B, Raoult D, Drancourt M. Archaea as emerging organisms in complex human microbiomes. Anaerobe. 2011; 17: 56–63. https://doi.org/10.1016/j.anaerobe.2011.03.001 PMID: 21420503

5. Carberry CA, Waters SM, Kenny DA, Creevey CJ. Rumen methanogenic genotypes differ in abundance according to host residual feed intake phenotype and diet type. Appl Env Microbiol. 2014; 80: 586–594.

6. Borrel G, McCann A, Deane J, Neto MC, Lynch DB, Brugère JF, et al. Genomics and metagenomics of trimethylamine-utilizing Archaea in the human gut microbiome. ISME J. 2017; 11: 2059–2074. https://doi.org/10.1038/ismej.2017.72 PMID: 28585938

7. Bang C, Weidenbach K, Gutsmann T, Heine H, Schmitz RA. The intestinal archaea Methanosphaera stadtmanae and Methanobrevibacter smithii activate human dendritic cells. PLoS One. 2014; 9: e99411. https://doi.org/10.1371/journal.pone.0099411 PMID: 24915454

8. Ghavami SB, Rostami E, Sephay AA, Shahrokh S, Balaii H, Aghdaei HA, et al. Alterations of the human gut Methanobrevibacter smithii as a biomarker for inflammatory bowel diseases. Microb Pathog. 2018; 117: 285–289. https://doi.org/10.1016/j.micpath.2018.01.029 PMID: 29477743

9. Mathur R, Barlow GM. Obesity and the microbiome. Expert Review of Gastroenterology and Hepatology. 2015. pp. 1087–1099. https://doi.org/10.1586/17474124.2015.1051029 PMID: 26082274

10. Hook SE, Wright A-DG, McBride BW. Methanogens: methane producers of the rumen and mitigation strategies. Archaea. 2010; 945785. https://doi.org/10.1155/2010/945785 PMID: 21253540

11. Poulsen M, Schwab C, Jensen BB, Engberg RM, Spang A, Canibe N, et al. Methylotrophic methanogenic Thermoplasmata implicated in reduced methane emissions from bovine rumen. Nat Commun. 2013; 4: 1428. https://doi.org/10.1038/ncomms2432 PMID: 23385573

12. Jarvis GN, Strömpl C, Burgess DM, Skillman LC, Moore ERB, Joblin KN. Isolation and identification of ruminal methanogens from grazing cattle. Curr Microbiol. 2000; 40: 327–332. https://doi.org/10.1007/s002849910065 PMID: 10706664

13. Friedman N, Jami E, Mizrahi I. Compositional and functional dynamics of the bovine rumen methanogenic community across different developmental stages. Environ Microbiol. 2017; 19: 3365–3373. https://doi.org/10.1111/1462-2920.13846 PMID: 28654196

14. Hansen EE, Lozupone CA, Rey FE, Wu M, Guruge JL, Narra A, et al. Pan-genome of the dominant human gut-associated archaeon, Methanobrevibacter smithii, studied in twins. Proc Natl Acad Sci. 2011; 108: 4599–4606. https://doi.org/10.1073/pnas.1000071108 PMID: 21317366

15. Thauer RK, Kaster AK, Seedorf H, Buckel W, Hedderich R. Methanogenic archaea: Ecologically relevant differences in energy conservation. Nat Rev Microbiol. 2008; 6: 579–591. https://doi.org/10.1038/nrmicro1931 PMID: 18587410

16. Morgavi DP, Forano E, Martin C, Newbold CJ. Microbial ecosystem and methanogenesis in ruminants. Animal. 2010; 4: 1024–1036. https://doi.org/10.1017/S1751731110000546 PMID: 22444607
17. Jackson BE, McInneray MJ. Anaerobic microbial metabolism can proceed close to thermodynamic limits. Nature. 2002; 415: 454–456. https://doi.org/10.1038/415454a PMID: 11807560

18. Gonzalez-Cabaleiro R, Lema JM, Rodriguez J, Kleerebezem R. Linking thermodynamics and kinetics to assess pathway reversibility in anaerobic bioprocesses. Energy Environ Sci. 2013; 6: 3780–3789.

19. Muñoz-Tamayo R, Laroche B, Walter E, Doré J, Leclerc M. Mathematical modelling of carbohydrate degradation by human colonic microbiota. J Theor Biol. 2010; 266: 189–201. https://doi.org/10.1016/j.jtbi.2010.05.040 PMID: 20561534

20. Van Wey AS, Lovatt SJ, Roy NC, Shorten PR. Determination of potential metabolic pathways of human intestinal bacteria by modeling growth kinetics from cross-feeding dynamics. Food Res Int. 2016; 88: 207–216.

21. Shoaire S, Ghaffari P, Kovatcheva-Datchary P, Mardinoglu A, Sen P, Pujos-Guillot E, et al. Quantifying Diet-Induced Metabolic Changes of the Human Gut Microbiome. Cell Metab. 2015; 22: 320–331. https://doi.org/10.1016/j.cmet.2015.07.001 PMID: 26244934

22. Kettle H, Louis P, Holtrop G, Duncan SH, Flint HJ. Modelling the emergent dynamics and major metabolites of the human colonic microbiota. Environ Microbiol. 2015; 17: 1615–1630. https://doi.org/10.1111/1462-2920.12599 PMID: 25142831

23. Heijnen JJ, Dijken JP. In search of a thermodynamic description of biomass yields for the chemotrophic growth of microorgansims. Biotechnol Bioengeneering. 1992; 39: 833–852.

24. Kleerebezem R, Van Loosdrecht MCM. A Generalized Method for Thermodynamic State Analysis of Environmental Systems. Crit Rev Environ Sci Technol. 2010; 40: 1–54.

25. Hoh CY, Cord-Ruwisch R. A practical kinetic model that considers endproduct inhibition in anaerobic digestion processes by including the equilibrium constant. Biotechnol Bioeng. 1996; 51: 597–604. https://doi.org/10.1002/(SICI)1097-0290(19960905)51:5<597::AID-BIT12>3.0.CO;2-F PMID: 18629824

26. Desmond-Le Quemener E, Bouchez T. A thermodynamic theory of microbial growth. Isme J. 2014; 8: 1747–1751. https://doi.org/10.1038/ismej.2014.7 PMID: 24522260

27. Großkopf T, Soyer OS. Microbial diversity arising from thermodynamic constraints. ISME J. 2016; 10: 2725–2733. https://doi.org/10.1038/ismej.2016.49 PMID: 27035705

28. Kohn R A, Boston RC. The Role of Thermodynamics in Controlling Rumen Metabolism. Model Nutr Util Farm Anim. 2000; 11–24.

29. Ungerfeld EM. A theoretical comparison between two ruminal electron sinks. Front Microbiol. 2013; 4.

30. Janssen PH. Influence of hydrogen on rumen methane formation and fermentation balances through microbial growth kinetics and fermentation thermodynamics. Anim Feed Sci Technol. 2010; 160: 1–22.

31. Van Lingen HJ, Plugge CM, Fadel JG, Kebreab E, Bannink A, Dijkstra J. Thermodynamic driving force of hydrogen on rumen microbial metabolism: A theoretical investigation. PLoS One. 2016; 11: e0161362. https://doi.org/10.1371/journal.pone.0161362 PMID: 27783615

32. Offner A, Sauvant D. Thermodynamic modeling of ruminal fermentations. 2006; 55: 343–365.

33. Ghimire S, Gregorini P, Hanigan MD. Evaluation of predictions of volatile fatty acid production rates by the Molly cow model. J Dairy Sci. 2014; 97: 354–362. https://doi.org/10.3168/jds.2012-6199 PMID: 24268399

34. Wolfe RS. Techniques for cultivating methanogens. Methods Enzymol. 2011; 494: 1–22. https://doi.org/10.1016/B978-0-12-385112-3.00001-9 PMID: 21402207

35. Ohene-Adjei S, Teather RM, Ivan M, Forster RJ. Postinoculation protozoan establishment and association patterns of methanogenic archaea in the ovine rumen. Appl Environ Microbiol. 2007; 73: 4609–4618. https://doi.org/10.1128/AEM.02687-06 PMID: 17513586

36. Popova M, Martin C, Eugene M, Mialon MM, Doreau M, Morgavi DP. Effect of fibre- and starch-rich finishing diets on methanogenic Archaea diversity and activity in the rumen of feedlot bulls. Anim Feed Sci Technol. 2011; 166–67: 113–121.

37. Bricheux G, Bonnet JL, Bohatier J, Morel JP, Morel-Desrosiers N. Microcalorimetry: a powerful and original tool for tracking the toxicity of a xenobiotic on Tetrahymena pyriformis. Ecotoxicol Env Saf. 2013; 98: 88–94.

38. Braissant O, Bonkat G, Wirz D, Bachmann A. Microbial growth and isothermal microcalorimetry: Growth models and their application to microcalorimetric data. Thermochim Acta. 2013; 555: 64–71.

39. Muñoz-Tamayo R, Giger-Reverdin S, Sauvant D. Mechanistic modelling of in vitro fermentation and methane production by rumen microbiota. Anim Feed Sci Technol. 2016; 220: 1–21.

40. Batstone DJ, Keller J, Angelidaki I, Kalyuzhnyi S V, Pavlostathis SG, Rozzi A, et al. Anaerobic Digestion Model No.1 (ADM1). IWA Task Group for Mathematical Modelling of Anaerobic Digestion Processes. IWA Publishing, London; 2002.
41. Pavlostathis SG, Miller TL, Wolin MJ. Cellulose Fermentation by Continuous Cultures of Ruminococcus-Albus and Methanobrevibacter-Smithii. Appl Microbiol Biotechnol. 1990; 33: 109-116.

42. Schauer NL, Ferry JG. Metabolism of formate in Methanobacterium formicicum. J Bacteriol. 1980; 142: 800-807. PMID: 6769911

43. Haydock AK, Porat I, Whitman WB, Leigh JA. Continuous culture of Methanococcus maripaludis under defined nutrient conditions. FEMS Microbiol Lett. 2004.

44. Walter E, Pronzato L. Identification of Parametric Models from Experimental Data. Springer, London; 1997.

45. Muñoz-Tamayo R, Puillet L, Daniel JB, Sauvant D, Martin O, Taghipoor M, et al. Review: To be or not to be an identifiable model. Is this a relevant question in animal science modelling? Animal. 2018; 12: 701-712. https://doi.org/10.1017/S1751731117002774 PMID: 29096725

46. Bellu G, Saccomani MP, Audoly S, D'Angio L. DAISY: A new software tool to test global identifiability of biological and physiological systems. Comput Methods Programs Biomed. 2007; 88: 52-61. https://doi.org/10.1016/j.cmpb.2007.07.002 PMID: 17707944

47. Vanrolleghem PA, Vandaele M, Dochain D. Practical identifiability of a biokinetic model of activated-sludge respiration. Water Res. 1995; 29: 2561-2570.

48. Muñoz-Tamayo R, Laroche B, Leclerc M, Walter E. IDEAS: A parameter identification toolbox with symbolic analysis of uncertainty and its application to biological modelling. IFAC Proceedings Volumes. 2009. pp. 1271-1276.

49. Lin LI. A concordance correlation-coefficient to evaluate reproducibility. Biometrics. 1989; 45: 255-268. PMID: 2720055

50. Schill NA, Liu JS, von Stockar U. Thermodynamic analysis of growth of Methanobacterium thermoautotrophicum. Biotechnol Bioeng. 1999; 64: 74-81. https://doi.org/10.1002/(sici)1097-0290(19990705)64:1<74::aid-bit8>3.0.co;2-3 PMID: 10397841

51. Liu JS, Marison IW, von Stockar U. Microbial growth by a net heat up-take: A calorimetric and thermodynamic study on acetotrophic methanogenesis by Methanosarcina barkeri. Biotechnol Bioeng. 2001; 75: 170-180. https://doi.org/10.1002/bit.1176 PMID: 11536139

52. Ruiz T, Bec A, Danger M, Koussoroplis A, Aguer J, Morel J, et al. A microcalorimetric approach for investigating stoichiometric constraints on the standard metabolic rate of a small invertebrate. Ecol Lett. 2018; 21: 1714-1722. https://doi.org/10.1111/ele.13137 PMID: 30151853

53. Janssen PH, Kirs M. Structure of the archaeal community of the rumen. Appl Environ Microbiol. 2008; 74: 3619-3625. https://doi.org/10.1128/AEM.02812-07 PMID: 18424540

54. von Stockar U, Larsson C, Marison IW. Calorimetry and energetic efficiencies in aerobic and anaerobic microbial growth. Pure Appl Chem. 1993; 65: 1889-1892.

55. Von Stockar U, Liu JS. Does microbial life always feed on negative entropy? Thermodynamic analysis of microbial growth. Biochim Biophys Acta-Bioenerg. 1999; 1412: 191-211.

56. Hardin G. The competitive exclusion principle. Science (80-). 1960; 131: 1292-1297.

57. Lynch TA, Wang Y, van Brunt B, Pacheco D, Janssen PH. Modelling thermodynamic feedback on the metabolism of hydrogenotrophic methanogens. J Theor Biol. 2019; 477: 14-23. https://doi.org/10.1016/j.jtbi.2019.05.018 PMID: 31150665

58. Udén P, Rounsaville TR, Wiggans GR, Van Soest PJ. The measurement of liquid and solid digesta retention in ruminants, equines and rabbits given timothy (Phleum pratense) hay. Br J Nutr. 1982; 48: 329-339. https://doi.org/10.1079/bjn19820117 PMID: 6810917

59. Ng F, Kittelmann S, Patchett ML, Attwood GT, Janssen PH, Rakonjac J, et al. An adhesin from hydrogen-utilizing rumen methanogen Methanobrevibacter ruminantium M1 binds a broad range of hydrogen-producing microorganisms. Env Microbiol. 2016; 18: 3010-3021.

60. Samuel BS, Hansen EE, Manchester JK, Coutinho PM, Henrissat B, Fulton R, et al. Genomic and metabolic adaptations of Methanobrevibacter smithii to the human gut. Proc Natl Acad Sci. 2007; 104: 10643-10648. https://doi.org/10.1073/pnas.0704189104 PMID: 17563350

61. Kelly WJ, Leahy SC, Li D, Perry R, Lambie SC, Attwood GT, et al. The complete genome sequence of the rumen methanogen Methanobacterium formicicum BRM9. Stand Genomic Sci. 2014; 9: 15. https://doi.org/10.1186/1944-3277-9-15 PMID: 25780506

62. Hutchinson GE. The paradox of the plankton. Am Nat. 1961; 95: 137-145.

63. Bernalier A, Lelait M, Rochet V, Grivet JP, Gibson GR, Durand M. Acetogenesis from H2 and CO2 by methane- and non-methane-producing human colonic bacterial communities. FEMS Microbiol Ecol. 1996; 19: 193-202.

64. Nava GM, Carbonero F, Croix JA, Greenberg E, Gaskins HR. Abundance and diversity of mucosa-associated hydrogenotrophic microbes in the healthy human colon. ISME J. 2012; 6: 57-70.
65. Flint HJ, Duncan SH, Scott KP. Interactions and competition within the microbial community of the human colon : links between diet and health. Environ Microbiol. 2007; 9: 1101-1111. https://doi.org/10.1111/j.1462-2920.2007.01281.x PMID: 17472627

66. Huws SA, Creevey CJ, Oyama LB, Mizrahi I, Denman SE, Popova M, et al. Addressing global ruminant agricultural challenges through understanding the rumen microbiome: past, present, and future. Front Microbiol. 2018; 9: 2161. https://doi.org/10.3389/fmicb.2018.02161 PMID: 30319557

67. Pfeiffer T, Schuster S, Bonhoeffer S. Cooperation and competition in the evolution of ATP-producing Pathways. Science (80-). 2001; 292: 50-507.

68. Vandermeer J, Evans MA, Foster P, Höök T, Reiskind M, Wund M. Increased competition may promote species coexistence. Proc Natl Acad Sci. 2002; 99: 8731-8736. https://doi.org/10.1073/pnas.142073599 PMID: 12070354

69. MacLean RC, Gudelj I. Resource competition and social conflict in experimental populations of yeast. Nature. 2006; 441: 498-501. https://doi.org/10.1038/nature04624 PMID: 16724064

70. Rapaport A, Dochain D, Harmand J. Long run coexistence in the chemostat with multiple species. J Theor Biol. 2009; 257: 252-259. https://doi.org/10.1016/j.jtbi.2008.11.015 PMID: 19111560

71. Grognard F, Masci P, Benoît E, Bernard O. Competition between phytoplankton and bacteria: exclusion and coexistence. J Math Biol. 2015; 70: 959-1006. https://doi.org/10.1007/s00285-014-0783-x PMID: 24748458

72. Hsu SB, Hubbell S, Waltman P. A mathematical theory for single-nutrient competition in continuous cultures of micro-organisms. SIAM J Appl Math. 1977; 32: 366-383.

73. Widder S, Allen RJ, Pfeiffer T, Curtis TP, Wiuf C, Sloan WT, et al. Challenges in microbial ecology: Building predictive understanding of community function and dynamics. ISME J. 2016; 10: 2557-2568. https://doi.org/10.1038/ismej.2016.45 PMID: 27022995

74. Muñoz-Tamayo R, Ramirez Agudelo JF, Dewhurst RJ, Miller G, Vernon T, Kettle H. A parsimonious software sensor for estimating the individual dynamic pattern of methane emissions from cattle. Animal. 2019; 13: 1180-1187. https://doi.org/10.1017/S1751731118002550 PMID: 30333069

75. Ou JZ, Yao CK, Rotbart A, Muir JG, Gibson PR, Kalantar-zadeh K. Human intestinal gas measurement systems: In vitro fermentation and gas capsules. Trends Biotechnol. 2015; 33: 208-213. https://doi.org/10.1016/j.tibtech.2015.02.002 PMID: 25772639