Universality of Thermodynamic Constants Governing Biological Growth Rates

Background Mathematical models exist that quantify the effect of temperature on poikilotherm growth rate. One family of such models assumes a single rate-limiting ‘master reaction’ using terms describing the temperature-dependent denaturation of the reaction's enzyme. We consider whether such a model can describe growth in each domain of life. Methodology/Principal Findings A new model based on this assumption and using a hierarchical Bayesian approach fits simultaneously 95 data sets for temperature-related growth rates of diverse microorganisms from all three domains of life, Bacteria, Archaea and Eukarya. Remarkably, the model produces credible estimates of fundamental thermodynamic parameters describing protein thermal stability predicted over 20 years ago. Conclusions/Significance The analysis lends support to the concept of universal thermodynamic limits to microbial growth rate dictated by protein thermal stability that in turn govern biological rates. This suggests that the thermal stability of proteins is a unifying property in the evolution and adaptation of life on earth. The fundamental nature of this conclusion has importance for many fields of study including microbiology, protein chemistry, thermal biology, and ecological theory including, for example, the influence of the vast microbial biomass and activity in the biosphere that is poorly described in current climate models.


Introduction
Temperature governs the rate of chemical reactions including those controlling the development and decline of life on earth from individual cells to complex populations. This has been so across geological time and will continue to be so as the world faces unprecedented uncertainty reflected in popular, political and scientific debate about the impact of temperature on the biosphere. Models to describe the effect of temperature on biological systems reliably are becoming increasingly important. One such family of models [1][2] assumes a single rate-limiting, enzyme-catalyzed 'master reaction' using an Arrhenius form modified by terms that describe the temperature-dependent denaturation of that enzyme. Earlier [3] we presented such a model that specifically describes the Gibbs free-energy change upon folding/unfolding of the putative master enzyme as a function of temperature and individually fitted it to 35 bacterial strains, but the approach required three parameters to be held fixed (T Ã H , T Ã S , and DS Ã ). Here we begin with the same assumption but use a Bayesian hierarchical modeling approach to analyze simultaneously 95 data sets from all three domains of life, Bacteria, Archaea and Eukarya [4], and, with only minimal constraints, estimate all thermodynamic parameters, some of which were previously only accessible by experiment. Importantly, the model implies that the temperature-related growth rate of all poikilotherm life behaves as if controlled by a single-enzyme system and the growth range is constrained by the denaturation of that enzyme.

Data
The data comprise standardized growth rates (or rates of metabolism in some cases) of 95 strains from 19 Bacteria, 19 Archaea, and 13 unicellular Eukarya species (Table S1). Here we use the yeasts to represent the domain Eukarya. The data were obtained by an ad hoc process, but are sufficiently diverse in origin and species to enable us to test the thermodynamic model's validity for physiologically different organisms and included observations at temperatures below and above the maximum growth rate. They include psychrophiles (e.g. Gelidibacter sp.), psychrotrophs (e.g. Shewanella gelidimarina), mesophiles (e.g. Escherichia coli), thermophiles (e.g. Streptococcus thermophilus), acidophiles (e.g. Ferroplasma acidiphilum), halophiles (e.g. Haloarcula vallismortis) and haloalkaliphiles (e.g. Natronococcus occultus).

Model Structure
Our model (equations 1-2) assumes that the growth rate is governed by a single enzymic rate-controlling reaction system that is limiting under all conditions. The quantity F is the growth rate given the temperature and the values of the parameters. The numerator is essentially an Arrhenius model that describes the rate of the putative enzyme-catalyzed rate-controlling reaction (RCR) as a function of temperature while the denominator (equation 2) models the change in expected rate due to the effects of temperature on the conformation and activity of the putative enzyme catalyzing that reaction and is essentially the reciprocal of the probability, as a function of temperature, that the enzyme is in a catalytically active conformation. where: In the model R is the gas constant (8.314 J/K mol); c is a scaling constant; DH { A the enthalpy of activation (J/mol) of the rate limiting reaction; T is the temperature (K); DC p is the heat capacity change (J/K mol-amino acid-residue) upon denaturation of the RCR; n is the number of amino acid residues in the RCR; DH Ã (J/mol-amino acid-residue) is the enthalpy change at T Ã H , the convergence temperature for enthalpy (K) of protein unfolding; and DS Ã (J/K mol-amino-acid-residue) is the entropy change at T Ã S , the convergence temperature for entropy (K) for protein unfolding. We derive two further quantities. One is the average number of non-polar hydrogen atoms per amino acid residue [5]: The DC p is associated with the energy required to reorganise the water molecules surrounding the protein [6] and increases with the non-polar accessible area of the molecule [7], as measured by N ch . The other is the temperature at which D in equation 2 reaches a minimum and at which denaturation is minimized [3]: DC p which provides an index of temperature adaptation of the organism. While the temperature of maximum growth could be used as an index, T mes is a more natural choice since it is the temperature at which the enzymic rate-controlling reaction system has evolved to be maximally active. At higher temperatures the growth rate may increase, but enzymic activity declines.
We allow four parameters to have values specific to each strain j~1,2, . . . 95 ð Þ : c j , DH { A j, n j , and DC p j . We assume the strain parameters to be Gaussian distributed with means specific to their domains: Bacteria, Archaea, or Eukarya. The remaining parameters (DH Ã , DS Ã , T Ã H , and T Ã S ) have values common to all strains and are thought not to depend on the individual biochemistry of each strain, but describe protein thermal stability limits [8][9][10]. In earlier analyses (data not shown) when data sets were grouped according to domain, estimates of the four thermodynamic constants related to protein stability were not significantly different between data sets and subsequently these were assumed to be universal, having the same values across all domains. This assumption is not disproven by the results and is supported by their narrow credible intervals. The common and domain parameters are each assigned a uniform prior with limits informed by the biochemistry literature. We assume that the square root of the observed growth rate has a Gaussian distribution with a mean given by the modeled value, F , and with an unknown precision (reciprocal variance).

Implementation
We use a Bayesian approach to allow for uncertainty in measurement and parameters to be incorporated in a natural way through the appropriate prior specification. We assign truncated normal priors to the strain parameters in which the means are specific to the domain of each, ð Þ is the domain (Bacteria, Archaea, Eukarya) for strain j. The exp t h ð Þ is the strain precision and models the variation between the strain parameters about the domain parameters, h d j ð Þ . The domain means and the t h are assigned uniform priors with limits informed by the biochemistry literature with the exception of c which is assigned a vague prior. The common parameters are each assigned a uniform prior with limits informed by the biochemistry literature. Finally, the observational precision is assigned a gamma distribution, t*C 0:001,0:001 ð Þ . Prior specifications are documented in Table S2. Inference is obtained in the form of posterior means and variances using Markov Chain Monte Carlo (MCMC) simulation [11]. We choose to update the parameters of each strain as a block using Haario updates [12]. We also use Haario updates [12] for each set of domain parameters and the strain parameter precisions. We use separate Metropolis-Hastings [13] updates for each of the common parameters but to ensure rapid mixing we also use a Metropolis-coupled MCMC approach [14] based on 60 chains. During each iteration we randomly select a pair of chains i,j ð Þ and swap their common parameters as a block with probability min 1, keeping other parameters unchanged. The model is run for 600,000 iterations and the last 50% of iterations are retained for further analysis. We summarize parameters using posterior block means, standard deviations, and 95% highest posterior density intervals (HPDI). A 95% HPDI is the shortest interval that contains a parameter with 95% probability. We determine the most likely ordering of the parameters in terms of their posterior probabilities by calculating the mean number of MCMC iterations for which each alternative arises.

Results
We obtain good fits for the growth curves of all strains as shown in Figure 1, and in more detail in Figures S1, S2, S3, S4, S5, S6. We summarize the strain parameters in Table S1 and Figure 2, common parameters in Table 1, and domain parameters in Table 2. Some strains appear outlying since they are below the 2.5 th or above the 97.  Figure 3 in which most of the strains with outlying parameters also have maxima that are either much lower than other strains (e.g. strains 20, 33, 35), or much higher (e.g. strains 48, 49, 68). This is also shown by the values of T mes for the strains 33, 35, 20, which are lower than other strains ( Figure 2). Some outliers may be explained when it is realized that those strains are at the extremes of thermal adaptation, such as strain 46, a thermophile, strain 48, a thermoacidophile, or strain 61, a psychrophile. In Figure 2, the strains 3, 19, 20, 33-35, appear to be outliers for T mes ; this is likely due to a lack of data in the lower temperature region, which also resulted in their wider credible intervals. Other variations in the data presented in Figure 2 are explicable after careful examination of Figures 1, 3, 4, and S1, S2, S3, S4, S5, S6.
On the right side axis of Figure  A , those parameters reflect the thermal stability of the RCR. The Eukarya strains included in this analysis are predominantly mesophilic, the Bacteria include more psychrophilic or psychrotrophic strains than the Archaea strains used, and the latter datasets include more thermophilic strains than the Bacteria, which included only two mildly thermophilic species, S. thermophilus and Acidimicrobium ferrooxidans, the former represented by multiple datasets (37)(38)(39)(40)(41)(42)(43)(44)(45)(46).
The remaining parameter, n, is largest in Eukarya, intermediate in Bacteria, and least in Archaea: P(E.B.A) = 0.98. Complete genome sequences give mean protein lengths of 449 for Eukarya, 330 for Bacteria, and 270 for Archaea [15]. With the most general assumptions our model is able to obtain the same ordering ( Table 2).
For those species represented by at least three datasets we calculated the percentage deviation for the strain parameters (Table S3): n (37.5%), DH { A (23.6%), N ch (5.7%), DC p (1.9%), T mes (0.5%). While the large variability for some parameters may be expected, the very low value for others such as DC p and T mes are of considerable interest. This can also be seen in Figure 2 where we use vertical shading to delimit species. The DC p and T mes appear almost species-specific, as is, to a lesser extent, N ch . The parameter T mes is mathematically related to DC p and so shares its low variability.
The  [24] and conceivably may relate specifically to yeast species used in wine production, rather than to the Eukarya as a whole.
The quantity P~1=D, the probability that the enzyme is in its native (catalytically active) state [3], is shown in Figure 4 for each strain. Some strains have narrower curves indicating a more limited range of temperatures within which the protein is in the native state, but others are ''flat-topped'' indicating a more extensive range. As can be seen from the extrapolated curves the optimal stability region is fully observed for most strains. The squat curves exhibited by the strains that do not reach maximum probability suggest that those strains were not grown under optimal conditions. As shown in Table 1, the 95% HPDI for T Ã H and T Ã S do not overlap, emphasizing that their values are not identical. Equivalently, the 95% HPDI for their difference, (216, 213), does not include zero. In fact, they are quite close to the values suggested by others [8][9][10], 373 K and 385 K, respectively, and which are contained within their 95% HPDI.

Whereas previous work [3] individually analyzed 35 Bacteria strains and held 3 parameters fixed (T Ã
H , T Ã S , and DS Ã ) we were, with only minimal constraints, able to analyze all 95 Archaea, Bacteria, and Eukarya strains to provide estimates of all strain, domain and common parameters simultaneously, and obtain reasonable estimates of thermodynamic parameters previously estimated decades ago by experiment [9][10]. The model successfully described the growth of unicellular organisms in all domains of life. We were able to do this by utilizing the MCMC methods to greatly simplify the computation involved with a Bayesian framework, and a hierarchical model that permits those parameters that are less well informed to 'borrow strength' from elsewhere in the model and so improve estimation.
While it is possible that the reaction system is more complex and only appears equivalent to a single rate-controlling system, single enzymes have been shown to be growth rate-controlling [25][26].  Table S1. Bacterial strains are shown as green circles, archaeal strains as blue squares, and eukaryote strains as red diamonds. doi:10.1371/journal.pone.0032003.g001 The model results lend support to the 'master reaction' model in which an enzyme-catalyzed reaction limits growth rate at all temperatures. The model indicates that, all other things being equal, thermodynamic properties of protein hydration govern biological rates and these are consistent in all forms of unicellular life. In turn this suggests that the thermal stability of proteins is a fundamental property in the evolution and adaptation of life on earth and has importance for many fields of study including microbiology, protein chemistry, thermal biology, and ecological theory, including, for example, the influence of the vast microbial biomass and activity in the biosphere that is poorly described in current climate models [27].  Since this is only the debut of this model we anticipate further developments and extensions. Our results appear consistent with the assumption that growth is controlled by a single rate-limiting, enzyme-catalyzed 'master reaction' using an Arrhenius form modified by terms that describe the temperature-dependent denaturation of that enzyme. While there may be other explanations for our results, candidates for a single-enzyme ratecontrolling mechanism common to all life can be conjectured. There may be processes that occur in parallel leading up to cell division, but it is also plausible that there is a single process dictated by a single enzyme (or other bio-catalyst) that could still be rate-limiting. There are various candidate biosynthetic processes common to all cells. Alternatively, anabolism or biosynthesis could be limited by the rate of energy generation e.g, via respiration. Regardless of whether there is a single ratelimiting, enzyme-catalyzed 'master reaction', or rate limiting process/pathway, there remains the possibility that some model parameters are more consequential at some temperatures than others. Whether the activation energy dominates the suboptimal range while denaturation is more important in the super optimal range is a topic for further investigation.
While we chose to adopt a model structure that recognises the three domains of life [4], the need for differentiation into three domain parameter sets remains to be explicitly considered. In the current study this was complicated by the availability of data. In particular, our data sets were not balanced in that they had  disproportionate representation of thermal adaptation groups (mesophiles, thermophiles, and so on). Further work is needed to examine whether an alternative, and perhaps simpler, domain parameter structure may be as effective. In addition, there is a need to consider if there are systematic differences of thermodynamic parameters between thermal adaptation groups, either within or between domains.
The results of our current model indicate that the activation energy of the putative 'master reaction' varies among strains and probably among species. This is in contradiction to the 'Metabolic Theory of Ecology' (MTE) [28][29] which describes metabolism in terms of an Arrhenius temperature dependence and assumes the same activation energy for all organisms. Other recent work also supports the idea that activation energy is not invariant [30][31] or interprets the variation of activation energy as being the result of protein denaturation [32], as in this study.
We only considered the temperature dependent growth of unicellular poikilothermic organisms, but it may be possible to apply the model directly to simple multicellular organisms or the growth of anatomical structures [33]. For larger organisms an allometric power function of body mass could be used, analogous to that also assumed in the MTE [28][29]34], whilst also remaining cognizant that the exponent in the function may not be the same for all forms of life [35]. While the current model well described temperature-dependent growth across the entire biokinetic range, practical application in environmental science is likely to require further extension to allow for additional factors such as nutrient and water availability and other environmental stressors. The extension to multicellularity and other environmental variables would broaden the model to a much wider ecological context. These would appear to be necessary steps before examining applications, e.g. in climate change, oceanography and soil carbon studies.
Many of these areas will have significant impact on aspirations of an improved human condition, especially for currently disadvantaged populations, balanced with the goal of maintaining biological diversity and activity in the long term. It has not escaped our attention that our work immediately suggests the possibility of temperature adaptation models applicable to wide ranging forms of life and life systems. Figure S1 Detailed fits for strains 1-16. Shown are the observed growth rate data using symbols and fitted curves for strains [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. Observed data are shown as green circles. All are strains of Bacteria. The fitted curves are calculated using the mean posterior parameter estimates and extend beyond the observed temperature range by 62.5u. (TIF)  Figure S2 Detailed fits for strains 17-32. Shown are the observed growth rate data using symbols and fitted curves for strain [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. Observed data are shown as green circles. All are strains of Bacteria. The fitted curves are calculated using the mean posterior parameter estimates and extend beyond the observed temperature range by 62.5u. (TIF) Figure S3 Detailed fits for strains [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]. Shown are the observed growth rate data using symbols and fitted curves for strain [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]. Observed data are shown as green circles for strains of Bacteria and blue squares for strains of Archaea. The fitted curves are calculated using the mean posterior parameter estimates and extend beyond the observed temperature range by 62.5u. (TIF) Figure S4 Detailed fits for strains 49-64. Shown are the observed growth rate data using symbols and fitted curves for strain 49-64. Observed data are shown as blue squares for strains of Archaea. The fitted curves are calculated using the mean posterior parameter estimates and extend beyond the observed temperature range by 62.5u. (TIF) Figure S5 Detailed fits for strains 65-80. Shown are the observed growth rate data using symbols and fitted curves for strain 65-80. Observed data are shown as blue squares for strains of Archaea and red diamonds for strains of Eukarya. The fitted curves are calculated using the mean posterior parameter estimates and extend beyond the observed temperature range by 62.5u. (TIF) Figure S6 Detailed fits for strains 81-95. Shown are the observed growth rate data using symbols and fitted curves for strain 81-95. Observed data are shown as red diamonds. All are strains of Eukarya. The fitted curves are calculated using the mean posterior parameter estimates and extend beyond the observed temperature range by 62.5u. (TIF)