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

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 CO2, 0.98 for H2 and 0.97 for CH4. 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.

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 2 , 0.98 for H 2 and 0.97 for CH 4 . 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. PLOS

Introduction
Methanogenic archaea inhabit the gastro-intestinal tract of mammals where they have established syntrophic interactions within the microbial community [1][2][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 Methanomassillicoccus 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). Cytochromelacking 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][20][21][22], although ttheoretical 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][26][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 thermodynamicbased 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 [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 ruminantium 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 computed 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 continuous and simultaneous recording as a function of time of the heat flux produced by 12 samples. The bath temperature was set at 39˚C; its long-term stability was better than ± 1x10 -4˚C over 24h. Each minicalorimeter was electrically calibrated. The specific disposable 4 mL microcalorimetric glass ampoules capped with butyl rubber stoppers and sealed with aluminium crimps were filled with 1.75 mL of Balch growth media and overpressed with 2.5 Pa of H 2 /CO 2 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 beginning of the microcalorimetric measurements. Actively growing cultures of methanogens (OD 660 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 inoculation. 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 microcalorimetric 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 -1 ). Analogously, the heat flux follows an exponential behaviour determined by the parameter μ c as described by [37,38].
The growth rate constant μ c can be determined by fitting the exponential part of the 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 2 /CO 2 mixture in the gas phase diffuses to the liquid phase. The H 2 and CO 2 in the liquid phase are further utilized by the pure culture to produce CH 4 . 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 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�(M � 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 where μ is the growth rate (h -1 ), μ 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 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 (ΔG r,c ) when the reaction is exergonic (Δ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 r,destock ) when the reaction is exergonic (ΔG r,destock <0) or zero otherwise. Finally, the dissipated exergy E dis 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 1 ) for methane production and one anabolic reaction (R 2 ) 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 5 H 7 O 2 N [40].
In the model, the stoichiometry of the reactions is taken into account via the parameters Y, Y CO 2 , Y CH 4 , which are the yield factors (mol/mol) of microbial biomass, CO 2 and CH 4 . 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 2 , using the methane-based molar growth yield of 2.8 g biomass/ mol CH 4 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 2 using the methane-based molar growth yield of 3.5 g biomass/mol CH 4 reported by Schauer and Ferry [42]. The fraction of H 2 utilized for microbial growth (reaction R 2 ) is defined by the yield factor Y (mol of microbial biomass/ mol of H 2 ). Now, let f be the fraction of H 2 used for the catabolic reaction R 1 . Reaction R 2 tells us that for every 10 moles of H 2 used in R 2 , we get 1 mol of microbial biomass. Hence, it follows that If Y is known, the fraction f can be computed from Eq (9). The yield factors of CO 2 and CH 4 can be expressed as functions of the the fraction f: The model has two physicochemical parameters (k L a, K H;CO 2 ) and four biological parameters (μ max , K s , Y, k d ). The initial condition for s CO 2 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 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 -1 ) is the dilution rate of the methanogens and b (h -1 ) 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 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 μ 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 μ max was equal to the specific rate constant μ 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 1 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 2 , CH 4 , CO 2 ). The Lin's concordance correlation coefficient (CCC) [49] was computed to quantify the agreement between the observations and model predictions.

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.  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 (ΔH m ), Gibbs energy (ΔG m ) and entropy (Δ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.66J. 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 μ max for each methanogen corresponds to   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.

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 2 . 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  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 (μ 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 2 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 μ 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 μ 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 2 /CO 2 is characterized by large heat production. The growth is highly exothermic, with a ΔH m value that largely exceeds the values found when other energy substrates are used. The enthalpy change ΔH m , which is more negative than the Gibbs energy change ΔG m , largely controls the process. Growth on H 2 /CO 2 is also characterized by a negative entropic contribution TΔ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ΔS m corresponds in fact to the balance between the final state and the initial state of the process, that is Methanogenesis on H 2 /CO 2 is particular because the final state of its catabolic reaction (1 mol CH 4 + 2 mol H 2 O) involves a smaller number of moles than the initial state (4 mol H 2 + 1 mol CO 2 ), which results in a significant loss of entropy during the process. For spontaneous growth in such a case, the ΔH m must not only contribute to the driving force but must also compensate the growth-unfavourable TΔS m , which means that ΔH m must be much more negative than ΔG m [55]. For this reason, methanogenesis on H 2 /CO 2 , 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 (ΔH m � 0 and TΔS m > 0), whereas fermentative metabolism is mainly entropy-driven (ΔH m < 0 and TΔS m � 0). Methanogenesis on H 2 /CO 2 is enthalpy-driven but entropy-retarded (ΔH m � 0 and TΔS m < 0), whereas methanogenesis on acetate is entropy-driven but enthalpy-retarded (ΔH m > 0 and TΔS m � 0). In the present case, the highly exothermic growth of M. ruminantium, M. smithii and M. formicium on H 2 /CO 2 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ΔS m and ΔH m 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 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 -1 (retention time = 48 h) and D = 0.04 h -1 (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] [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 -1 and the hydrogen flux production q H 2 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 -1 ). 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 coexistence 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, methanogenic diversity in the rumen results from multiple factors that include pH sensitivity, the association 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 2 via different pathways [63][64][65]. Both in the human gut and in the rumen, microbes grow in association with biofilms that form a polymerbased 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][68][69][70][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