HetMM: A Michaelis-Menten model for non-homogeneous enzyme mixtures

Summary The Michaelis-Menten model requires its reaction velocities to come from a preparation of homogeneous enzymes, with identical or near-identical catalytic activities. However, this condition is not always met. We introduce a kinetic model that relaxes this requirement, by assuming there are an unknown number of enzyme species drawn from a probability distribution whose standard deviation is estimated. Through simulation studies, we demonstrate the method accurately discriminates between homogeneous and heterogeneous data, even with moderate levels of experimental error. We applied this model to three homogeneous and three heterogeneous biological systems, showing that the standard and heterogeneous models outperform respectively. Lastly, we show that heterogeneity is not readily distinguished from negatively cooperative binding under the Hill model. These two distinct attributes—inequality in catalytic ability and interference between binding sites—yield similar Michaelis-Menten curves that are not readily resolved without further experimentation. Our user-friendly software package allows homogeneity testing and parameter estimation.


INTRODUCTION
One hundred and ten years after its initial publication, 1 the equation devised by Leonor Michaelis and Maud Menten is still routinely applied to predicting and explaining biochemical systems, to the point where it is almost synonymous with enzyme kinetics.And like all models, the Michaelis-Menten (MM) model makes several assumptions about the conditions that generated the observed data (see review by Srinivasan 2022 2 ).It is assumed that the enzyme concentration is many orders of magnitude less than its substrate, 3 and the latter is assumed to remain stationary for the observed duration of the reaction, 4 as does the concentration of the enzyme-substrate complex in the quasi-steady state that characterizes the model. 5nother fundamental assumption, this one so intuitively obvious that it can be easily overlooked, is that the catalysts (and reactants) are homogeneous.That is, the solution is composed of an ensemble of catalysts, each of which displays identical, or near-identical, catalytic activity.However, this condition is often violated in biological systems, even in purified enzyme preparations. 6,7There may exist multiple enzyme species in the preparation, each with distinct catalytic rates k cat and dissociation constants K D .In such a scenario, it would be preferable to characterize the kinetics of the catalyst, not as a single entity, but as a mixture of catalytic agents, whose catalytic properties are described by some probability distribution (Figure 1).Indeed, there are numerous cases of nonhomogeneous biocatalytic systems, both in vitro and in vivo.Many genomes express isozymes; enzymes with distinct sequences, but common ancestry and common catalytic activity, such as the multiple aminoacyl-tRNA synthetase gene duplicates that coexist in bacterial cells. 8In some cases, purified enzyme preparations contain multiple isozymes, such as cytochrome oxidase. 7,9But even a single gene can give rise to a heterogeneous mixture of protein products, through transcriptional or translational errors, alternative initiation or splicing, or post-translational modification, among other processes.Otherwise-rare occurrences can be amplified by the environment, for example thermophilic Archaea can adapt to lower temperatures by increasing the rate of methionine-for-leucine mistranslations. 10Primordial proteins in particular were most likely produced by a low-fidelity and ambiguous genetic code, and hence would have been highly variable. 11,12And even proteins with identical sequences, devoid of post-translational modifications, can still fold into a population of distinct structures 13,14 with a population of catalytic abilities, 15 or lack thereof -for instance ion channels 16 and molten globular enzymes. 17Any of these sources of variability could be further compounded by various degradation effects that may emerge from natural biological processes or from an ill-prepared in vitro sample.Catalase, for instance, is known as a ''suicide'' catalyst because it is degraded by its own substrate. 18That substrate is hydrogen peroxide, a reactive agent that oxidizes proteins; cysteine side chains in particular. 19,20urther examples of heterogeneity are illustrated by Brown et al. 2014. 7 On a coarser scale, heterogeneity may be an intrinsic property of the system that cannot be reduced to some mere idiosyncrasy of a particular gene or protein. 21This may be the case when studying cells, organelles, or multicellular organisms.When describing such a system with the MM model, the 'enzyme' would be a series of complex in vivo processes, and the 'enzyme preparation' may be living tissue.For instance, Devaux et al. 2023 22 applied the MM model to quantify oxygen consumption in fish brains.
Generally, heterogeneity is something to avoid in experimental design, but in some cases it may be desirable or even unavoidable.Suppose that not just one, but rather a large sample of enzymes, were to be co-expressed, co-purified, and co-assayed as part of a combinatorial expression library.The assay would measure the kinetic properties of the population of forms rather than just one enzyme species.For instance, experimental investigations into ancestral models typically operate on one enzyme sequence at a time, [23][24][25][26] but this approach could benefit from co-assaying a sample of reconstructions (such as a phylogenetic posterior distribution) to improve the robustness and statistical rigor of the results in a cost-affordable manner.
Some enzymes, such as aspartate kinase 27 and UDP-GlcNAc 2-epimerase, 28 display negative cooperativity between binding sites, where binding at one site interferes with the binding affinities at other sites.0][31] However, as discussed by Abeliovich 2005, negative cooperativity is not readily distinguished from a mixture of independent but variable binding sites. 32They propose the 1/N rule, where a Hill coefficient less than the reciprocal of the number of binding sites likely arose from heterogeneity, rather than negative cooperativity.
Clearly, there are many instances where the standard homogeneous model is unjustifiable, and many more instances where one may wish to test the validity of this assumption.In this article we describe an extension to the Michaelis-Menten model, in which it is assumed that there are an unknown number of catalytic species whose properties are independently drawn from an unknown probability distribution, and thus the enzyme preparation is heterogeneous.We provide a Bayesian model averaging framework, to test whether the standard assumption of homogeneity is adequate, or whether the heterogeneous model is preferred.Lastly, we explore the interactions between the heterogeneous and Hill models, and discuss the limitations of either approach.These methods are implemented in the user-friendly, open-source heterogeneous Michaelis-Menten (HetMM) package, which employs the BEAST 2 Bayesian inference engine. 33HetMM is available online at https://github.com/jordandouglas/HetMM.

Homogeneous and heterogeneous models
Consider the following reaction scheme between substrate A, enzyme E, and product P.

EA/
kcat E + P (Equation 1)  2) for substrate concentration a, where v is shorthand for V max (the maximum reaction rate) and k for the Michaelis constant K m (the substrate concentration where b r is half of V max ).Under the steady-state approximation, k is defined as: 3) The interpretation of k may change if any of the standard Michaelis-Menten assumptions 2 are not satisfied.Here, we will denote the natural logarithms of v and k by their uppercase symbols V = log v and K = log k.The transformed parameters V and K are not necessarily positive.Now, let us assume that V and K are independently sampled from normal distributions, or equivalently, v and k are independently sampled from log-normal distributions: 4) 5) Here, m ˛R and s > 0 are the respective means and standard deviations of V and K.These two normal distributions have probability densities f V and f K .Now let us marginalise across all values of K and V according to this distribution: b r ðaÞ = 8) for some function gðVÞ.As this function is independent of a, it is clear that any degree of variation in V cannot be resolved here without further information, such as the relative abundances of the species.Going forward, we will define v 0 as the mean V max across the species, weighted according to their unknown abundances.However, this simplification is not possible in the case of K.
dK (Equation 9) A numerical approximation of this integral is used.This new parameter, s K , describes the degree of variability in k.As s K approaches zero, the homogeneous and heterogeneous models become indistinguishable, and as s K grows beyond 1-2, the Michaelis constants k span several orders of magnitude (right panel of Figure 2).As confirmed in Figure 2, the Michaelis-Menten curve is sensitive to changes in s K when s K > 1. Whereas, when s K > 1, the assumption of homogeneity is often adequate, even if untrue.Overall, this model of heterogeneity assumes there are an unknown number of enzymes whose Michaelis constants fall onto a continuous spectrum of binding affinities and catalytic rates.In practice, the various k values of the mixture need not be log-normally distributed, however this flexible distribution is significantly more relaxed than the assumption of all catalytic species having identical k.This continuous model differs from the heterogeneity model by Brown et al. 2014, 7 which assumes a small number of individually homogeneous species (e.g., two species) that are weighted according to their relative abundances.

How heterogeneous is heterogeneous?
Before the distinction between homogeneity and heterogeneity makes any sense, we must first quantify how disparate the catalytic species should be so that their population is observably nonhomogeneous.Suppose, for instance, that an ensemble of enzymes were composed of several species, whose K m ranged from 1 to 1.1 units (if the species were to be characterized in isolation).In this case, the system would be effectively homogeneous.As demonstrated by Brown et al. 2014, the assumption of homogeneity is likely sufficient for the case of two species, provided that the ratio between their K m is no more than 20. 7Moreover, as shown in Figure 2 (bottom right), the standard deviation s K should be > 1 in order for the relative difference between the upper and lower extremes of k (among the population of enzymes) to be multiple orders of magnitude apart, and for the two models to have distinguishable properties.This logic forms the basis of our prior distribution for s K , which is shown in Table 1 along with the other priors.

Bayesian inference and model averaging
Parameter estimation is performed using Bayesian inference, allowing for robust estimation of parameters and their credible intervals, as well as the model indicator I m ˛f0; 1g, which governs the use of the homogeneous or heterogeneous MM model.This approach is known as Bayesian model averaging. 34,35It enables robust and probabilistic estimation of the model as if it were a parameter.Much like the Akaike and Bayesian information criteria, the method penalizes overparameterized models.Let D be the data measured from a kinetic assay, consisting of a vector of substrate concentrations a = ða 1 ; a 2 ; .; a n Þ and empirical reaction rates r = ðr 1 ;r 2 ;.;r n Þ.Given model parameters q, the likelihood is calculated by comparing observations r with their expected values b r , conditional on substrate concentration a, using the equation 10) where the experimental errors between expected b r ðaÞ and observed r values are assumed to follow a normal distribution, in log-space.This normal distribution, with probability density f e , has mean 0 and standard deviation e. Modeling experimental error in log-space accounts for the fact that reaction rates are non-negative, while also accommodating for the common scenario where the magnitude of error increases with substrate concentration, i.e., heteroscedasticity.Under this parameterization, e describes relative error and is therefore agnostic about measurement units.An error of e = 0:2, for example, has the same meaning regardless of whether velocities are in units of M per second or mM per second.This model consists of five parameters q = ðv 0 ;k 0 ;s K ;e;I m Þ.If I m = 0, then the enzymes are assumed to be homogeneous, in which case b r is calculated using Equation 2, where v = v 0 and k = k 0 .s K is not being used in the model.Whereas, if I m = 1, then the enzymes are assumed to be heterogeneous, in which case b r is calculated using Equation 9, where m K = log k 0 À 1 2 s 2 K , so that m K reflects the mean of k in real-space, rather than log-space.In this model, all five model parameters in q are being used.Lastly, the posterior density is pðqjDÞ f pðDjqÞ 3 pðqÞ: (Equation 11) Bayesian inference was run on these datasets, with both models (homogeneous and heterogeneous).In the second two cases (s K = 2; 3), the heterogeneous model is best, because pðI m = 1Þ is large.However, when s K = 1, the homogeneous and heterogeneous model could not be discriminated between, and therefore the model averaging favored the simpler model I m = 0.
Our prior distributions pðqÞ are summarised in Table 1.The posterior distribution is sampled using Markov chain Monte Carlo (MCMC), as implemented in BEAST 2. 33 Although BEAST 2 was designed for phylogenetic inference, its use of efficient proposal kernels, including Bactrian, 36 adaptable variance multivariate normal, 37 and adaptable operator sampler 38 kernels, makes it an attractive engine for the purposes of this model.MCMC chains are run until the effective sample sizes of all parameters exceed 200, as diagnosed by Tracer. 39astly, we wish to address a point of potential confusion concerning the varying use of log-normal distributions.We have used log-normal distributions in three different contexts here.First, in Equation 5, we assume that the properties of enzymes in the population are distributed in a log-normal fashion, i.e., log-normal ðm K ;s K Þ.Second, in Equation 10, we assume the experimental error is drawn from log-normal (0;e) -this distribution describes random error in the experimentation process.Lastly, in Table 1, we use the log-normal as prior distributions -these distributions describe information about any a priori expectations of a parameter, before performing any analyses.As shown in Figure 2, log-normal distributions are flexible, they have positive domains, and their shapes are determined by only one parameter, s, making them readily interpretable.

Hill model extension
We describe a heterogeneous extension to the Hill model.0][31] The nature and degree of this interaction is described by the Hill coefficient h > 0: b r ðaÞ = va h k h +a h (Equation 12) When h > 1, this indicates positive cooperativity between binding sites (giving a sigmoidal MM curve), and h > 1 indicates negative cooperativity (giving an MM-like curve with a sharpened hyperbola).When h = 1, there is no interaction, and the Hill model is equivalent to MM.Many have argued that the Hill model is merely a descriptive model, and lacks a strong theoretical basis. 31In any case, like the MM model, the Hill model assumes homogeneity.Thus, under the heterogeneous Hill model: b r ðaÞ = Z K va h e Kh +a h f K ðKÞ dK (Equation 13) The posterior distribution of this model q h is estimated using MCMC from the posterior density (Equation 11).However, there are two additional parameters in q h that are absent from the heterogeneous MM model.These two parameters govern whether the Hill coefficient Note that a Log-normal distribution's m and s are the mean and standard deviation in log-space.The mean in regular space is log m À s 2 2 .In a real-world example, the priors for v 0 and k 0 should be adjusted to reflect the units of measurement and other a priori information about the system being studied.
is positive, negative, or neutral.h is equal to h 0 when I h = À 1, 1 when I h = 0, and 1 h0 when I h = 1, where h 0 ˛ð0; 1Þ is the Hill coefficient in the case of negative cooperativity, and its reciprocal for positive cooperativity.This parameterization ensures symmetry in their prior distributions, and places the probability mass away from h = 1 so that the three models are distinct.Moreover, their prior distributions should reflect prior knowledge about the number of binding sites, and place the majority (or all) of their probability densities within the range À 1 N ; 1

Á
and ð1; NÞ, respectively, where N is the number of binding sites in each enzyme.Inferring the value of the model indicator I h fÀ 1; 0; 1g during MCMC can be used to test whether the Hill model is applicable to a given dataset, and if so, then the direction in which cooperativity acts.This is also known as Bayesian modeling averaging.

Validation on simulated data: The model is well-calibrated
In order to test how accurately our method can recover parameters from datasets generated by the same model, we simulated two hundred datasets and recovered the parameters on each of them using Bayesian MCMC.These datasets were simulated using parameters that were randomly sampled from the prior distribution.This experiment was designed to test (1) whether the true parameters can be recovered, and (2) whether the true model (homogeneous or heterogeneous) can be recovered.It would be undesirable if, for instance, the heterogeneous model was disproportionately selected even on homogeneous data.The conditions of these simulations mirrored a typical experimental setup, with three sets of reaction rates independently measured per substrate concentration, across seven different concentrations.These seven concentrations flanked the Michaelis constant k, which was sampled from the prior with a mean value of 3 units.These concentrations were 0.5, 2.5, 5, 10, 50, 250, and 1000 units.These results provide confidence that the method is able to recover the true parameters and model.In the case of continuous parameters, the true value was in the 95% credible interval approximately 95% of the time (see Figure 3 coverage).Moreover, the relative differences between lower and upper estimates were sufficiently small that their parameter estimates remained fairly informative (see Figure 3 span).Notably, V max and e had spans less than two-fold, while the Michaelis mean k 0 and standard deviation s K had the largest spans (43 and 3 3 ), suggesting they were more challenging to estimate.The spans are likely to widen with less available data.The true model indicator I m was identified, with over 95% support, the majority of the time.In some cases, the model could not be confidently resolved, where 0:05 < pðI m = 1Þ < 0:95.These simulation studies provided no evidence of systematic ''overfitting'', where the model with a larger dimensionality (i.e., the heterogeneous model) was favored because it captures random noise.Moreover, discrimination between the two models can be achieved even on datasets with moderate degrees of experimental error, where e(0:5.The model is therefore quite robust to random error -a standard deviation of e = 0:5 means that 95% of the measured reaction velocities at a given Two hundred datasets were simulated, and the parameters used to generate them were recovered with Bayesian MCMC.The coverage is the percentage of simulations where the true parameter lies in its 95% credible interval (blue lines), compared to the times when its estimate is wrong (black lines).The span is the relative difference between the upper and lower 95% interval points, geometrically-averaged across all replicates and rounded to 2 significant figures.These results show that coverage is close to 95% and that the spans are within an order of magnitude, thus providing confidence in the ability of the method to recover parameters.s K is not part of the homogeneous model, and therefore its coverage is only calculated when pðI m = 1Þ > 0:95.In the bottom middle panel, the estimated model is indicated if it has more than 0.95 posterior support, while '?' denotes uncertainty (i.e., 0:05 < pðI = 1Þ < 0:95).concentration can span up to 10-fold, and yet the true model (homogeneous/heterogeneous) can often still be recovered.We performed the same experiment for the homogeneous Hill model and found similar results (Figure 4).Namely, the two models (Hill and not-Hill) could be accurately discriminated between, the parameters could be recovered, and the 95% credible intervals were small enough to be useful.Moreover, detection of Hill binding was even more robust to experimental error than heterogeneity.][43] We also explored the possibility of resolving the newly introduced heterogeneous model from the standard Hill model.To do this, we simulated and performed Bayesian inference on a further 400 datasets, but this time they were simulated under 233 = 6 different models: homogeneous and heterogeneous, crossed with positive, negative, and neutral Hill.Half of these datasets were simulated at 7 substrate concentrations, and the other half with 14 concentrations (3 replicates per concentration).These results suggest that the negative Hill model, which is often interpreted as describing interference between binding sites, is in many cases, indistinguishable from the heterogeneous model, i.e., the two models are non-identifiable (Figure 5).Doubling the number of observations improved the resolution but did not solve the issue.In simple terms, a nonhomogeneous enzyme preparation may give the appearance of an interaction between binding sites, and a system whose binding sites are interacting may give the appearance of an inequality in their catalytic abilities.These two distinct mechanisms are not readily disentangled without further information, or a large volume of high-quality kinetic data.

Negative control: Homogeneous biological systems
We tested our method on three kinetic datasets that were obtained from putatively homogeneous enzyme preparations.In each case, we tested for both homogeneity and use of the Hill coefficient in a joint Bayesian analysis.First, we considered the original experiments carried out by Michaelis and Menten in 1913, 1 which have since been translated into English. 44They measured the hydrolysis of sucrose (into glucose and fructose) by an enzyme then-known as invertase (enzyme commission number EC. 3.2.1.26),across seven sucrose concentrations.Second, we considered the hydrolysis of o-nitrophenyl-b-D-galactopyranoside (ONPG) by b-galactosidase (EC.3.2.1.23).These data are used in the renz package for R, 45,46 describing eight experimental replicates, each one carried out by a different group of second-year college students.Third, we considered the human hemoglobin (Hb) dataset from Severinghaus 1979, 47 which has been described as a ''gold standard'' for cooperative binding data. 48The Hill model is often applied to modeling Hb, as it describes the interaction between its binding sites.These three datasets come with varying standards of experimental precision, as quantified by e (Table 2).The Michaelis-Menten dataset had the lowest experimental error (e = 0:033), while unsurprisingly, ONPG, which was aggregated from eight student groups, was comparatively noisy (e = 0:35).The 95% credible intervals were reasonably small on all three datasets-usually spanning less than two-fold-making for informative parameter estimates.In all three cases, the homogeneous model was selected (Table 2; Figure 6).Specifically, the ONPG and Hb datasets strongly rejected the heterogeneous model, with pðI m = 1Þ < 0:05, while MM was less confident.The Hill model was rejected by MM and ONPG, and selected for Hb, with h = 1:72, consistent with positive cooperativity between binding sites.These experiments further corroborate our simulation studies, confirming that our model, and its extra parameter, do not have a tendency to overfit to random noise in the data.

Positive control: Heterogeneous biological systems
Finally, we tested our method on three biological systems, each suspected to be nonhomogeneous (Figure 7).These datasets were tested for both homogeneity and use of the Hill coefficient.Where appropriate, the Hill coefficient was heavily constrained a priori within À 1 N ; N Á for N binding sites, in accordance with the 1/N rule. 32The first two biological experiments were performed under conditions where the enzyme concentration was significantly less than the substrate, and thus the quasi-steady state assumption of the MM model should be satisfied.
First, we considered an enzyme that is degraded by its own substrate.This enzyme is manganese peroxidase (MnP; EC. 1.11.1.13),and it catalyzes the oxidation of Mn 2+ to Mn 3+ in the presence of H 2 O 2 .However, H 2 O 2 is a reactive oxygen species, which decomposes into an even more reactive agent HO À in the presence of metal ions, leading to protein degradation. 19 Palma et al. 1997 showed that Mn 2+ exerted a positive effect on enzymatic activity, whereas excess levels of H 2 O 2 in fact lowered the reaction rate. 49We applied our method to initial oxidation rates at varying levels of Mn 2+ , at time zero.These reactions were initiated at four levels of H 2 O 2. 49 Our method rejected the assumption of homogeneity in this system.Interestingly, the standard deviation of Michaelis constants declined with increasing doses of H 2 O 2 (from s K = 4:43 for 100 mM down to 3.02 for 800 mM), and so did the probability of heterogeneity (from pðI m = 1Þ = 1 for 100 mM down to 0.92 for 800 mM).In the last case (800 mM), the heterogeneous model was preferred, but the homogeneous model was still adequate.These results are consistent with reactive products of hydrogen peroxide attacking the protein ensemble, leading to subpopulations of enzymes with severely impaired catalytic activities, while still remaining functional.Greater doses of H 2 O 2 provided a more uniform degradation effect.As MnP is monomeric with a single active site, 50 the Hill model was not considered.Variability in Michaelis constant is the best explanation for the non-MM behavior in this system.
Second, we considered a transmembrane ion channel that adopts ''open'' and ''closed'' conformations.This channel is the human cystic fibrosis transmembrane conductance regulator (CFTR) protein (EC.5.6.1.6),and it binds ATP in order to open the channel, enabling the flow of chloride-ions along their electrochemical gradient. 51,52The channel closes upon ATP hydrolysis.The catalytically active (open) state is promoted by phosphorylation of the R domain. 16As this is a complex multi-state system, it is reasonable to hypothesize that the active sites may take on a broad distribution of catalytic affinities, and therefore behave heterogeneously.CFTR contains two ATP binding sites, 53 and therefore its Hill coefficient should be no less than 0.5 and no more than 2 in order to explain cooperative binding.This constraint was reflected in our prior distributions.Indeed, our method favored the heterogeneous model, with h = 1, in both datasets by Liu et al.  2017 (phosphorylated and dephosphorylated CFTR).The Hill model was not selected and therefore heterogeneity is a much better Figure 5. Exploring the interaction between heterogeneous and Hill models Datasets were simulated under six models: homogeneous Michaelis-Menten (MM), homogeneous negative Hill (h > 1), homogeneous positive Hill (h > 1), heterogeneous Michaelis-Menten (Het), heterogeneous negative Hill (Het h > 1), and heterogeneous positive Hill (Het h > 1).Bayesian inference was performed on each dataset to estimate the model and its probability.Where any model was identified with greater than 50% posterior support, it is indicated accordingly on the y axis, and where no single model could be identified, it is labeled with a '?'.These datasets were simulated with three replicates of seven (left) and fourteen (right) concentrations of substrate, using parameters sampled from the prior.The labeled probabilities are conditional on the true model (i.e., the columns sum to 1).These results suggest that the negative Hill model and heterogeneous model are often non-identifiable and the correct model cannot always be recovered, even on larger datasets, corroborating the findings of Abeliovich 2005. 32able 2. Testing biological datasets for homogeneity and their use of the Hill coefficient

Dataset
V Shown are median parameter estimates (and 95% credible intervals), rounded to 3 significant figures.The estimates for s K and h 0 are omitted when the heterogeneous and Hill models are not being used, respectively.
explanation for the non-MM behavior displayed by CFTR than cooperative binding.Without invoking a multistate model of CFTR domain rearrangement, 16 this coarse-grained description of heterogeneity remains the best explanation.Third, we considered an in vivo assay on fish brains.Devaux et al. 2023 evaluated oxygen consumption in the brains of four species of triplefin fish. 22The efficiency of the electron transfer chain is likely to vary from mitochondrion-to-mitochondrion and cell-to-cell, and therefore a living tissue like this is likely to be intrinsically nonhomogeneous.In contrast to our previous systems, constraining the Hill coefficient to reflect the number of binding sites is no longer straightforward as this is a living tissue and cooperativity could come in many forms emerging from complex metabolic and regulatory pathways.Moreover, the quasi-steady state assumption of the MM model is often violated in living systems. 3Across the four fish brains, our method selected either the heterogeneous model, the Hill model with negative cooperativity, or both.Thus, it is unclear which explanation -heterogeneity, negative cooperativity, or the multi-step reaction chain -is most suitable without further information.
In all three nonhomogeneous systems, our Bayesian approach was unable to reliably estimate the Michaelis means k 0 , whose credible intervals spanned anywhere between 6-fold and over a thousand-fold, thereby rendering the estimates almost meaningless.In contrast, the remaining parameters were estimated more informatively, usually spanning less than three-fold.Our simulation studies did not observe any estimates spanning more than an order of magnitude, suggesting that further work is required to recover Michaelis constants from these systems, for example by developing new models or incorporating more data into the analysis.Estimating the Michaelis constants of heterogeneous enzymes remains an open challenge.

DISCUSSION
In this article we presented a model that tests for and captures heterogeneity inherent in the Michaelis constants (K m ) of enzymic systems.The method was validated using simulated data (Figure 3) and then applied to three biological datasets where it greatly outperformed the standard Michaelis-Menten model (Figure 6).This model features an additional parameter s K that describes the standard deviation of Michaelis constants between active sites, which are assumed to be independently drawn from a log-normal distribution.But despite the additional dimensionality, there is no evidence of this model systematically overfitting to random noise in the data.Moreover, the method is quite tolerant to random experimental error.The true model can often be recovered even when the observed reaction velocities span up to one order of magnitude across replicates, but of course less error is always preferred (Figure 3).Corroborating the results of Brown et al., 7 we showed that the standard assumption of homogeneity imposed by the Michaelis-Menten model is likely to be sufficient in most cases, and failing only under extreme conditions, such as when the different enzymes' Michaelis constants span orders of magnitude (Figure 2).Under these conditions, some enzymes must be significantly less proficient than others, but still active, otherwise their catalytic activities would not be detected.
Heterogeneity is often indistinguishable from negative-cooperation under the Hill model 32 ; the latter has been widely criticized for lacking a strong theoretical basis. 31There exists a twilight zone, in which independently heterogeneous enzymes give the appearance of being nonindependently interfering, and vice versa, and the two fundamentally distinct properties give rise to similar hyper-rectangular Michaelis-Menten curves.The two properties are not readily disentangled without further information, such as the number of binding sites per enzyme. 32This represents a fundamental limitation in both models, as well as the broader approach of inferring the behavior of molecular processes from two dimensional hyperbolic curves.
We tested six biological systems for their homogeneities and cooperative binding.These systems were (1) sucrose hydrolysis from Michaelis and Menten 1913, 1 (2) ONPG hydrolysis from second-year college students, 45 (3) oxygen binding by hemoglobin, 47 (4) an enzyme (MnP) that is degraded by its own substrate, 49 (5) an ion channel (CFTR) that adopts multiple conformations, 16 and (6) oxygen consumption in a living tissue (fish brain). 22The first three systems were anticipated to be homogeneous, and our method confirmed this.Whereas, the final three systems were suspected to be heterogeneous.Our results suggested that the MnP and CFTR systems were indeed heterogeneous, and their non-MM-like curves could not be explained by cooperativity between binding sites.Middle: membrane protein CFTR has heightened ATPase activity when open and phosphorylated, however its varying conformations and varying degrees of phosphorylation may perform this reaction at different effective rates. 16Bottom: the fish brain mitochondria, and all their various transport, regulatory, and metabolic processes, consume oxygen and produce reactive oxygen species (ROS) in a non-MM manner. 22e demonstrated how the assumption of homogeneity can be relaxed in the case of the Michaelis-Menten and Hill models.However, the approach is flexible and may readily be incorporated into other kinetic models, such as the total quasi-steady state approximation. 3,54lthough our model of heterogeneity has explanatory power, the approximation of heterogeneity is not based on any first principles (aside from the general observation that many things in nature are approximately normally distributed).In many cases, a specialized model would be preferred, such as a multi-state model for ion channels, 16 gene regulation, 55 or transcription elongation. 56Such models may also help in yielding informative estimates of the mean Michaelis constant, which could not be reliably recovered in our three biological heterogeneous systems.
Historically, scientists have taken great strides to actively avoid nonhomogeneous enzyme preparations, so as to eliminate as much background noise as possible.However, in some cases heterogeneous preparations may be preferable.For example, one may wish to purify a large sample of computationally modeled enzyme sequences, such as ancestral reconstructions, and screen them simultaneously for catalytic activity.Moreover, capturing variability may be essential for characterizing primordial proteins, which may have existed as populations of ''quasi-species'' 57 produced by an ambiguous genetic code. 11Our results here have paved the way for these kinds of experiments.Heterogeneity is everywhere in nature and identifying it is only the first challenge.The greater challenge lies in harnessing heterogeneity experimentally and determining when it is biologically meaningful. 21

Limitations of the study
This model has limitations.First, heterogeneity cannot be readily distinguished from negatively cooperative binding under the Hill model.Second, our assumption of heterogeneity is captured by a log-normal distribution, however, this approximation of the spread of K m is not based on any first principles.Third, K m is often non-identifiable in heterogeneous systems.Fourth, although the method captures variability in K m , it does not capture variability in V max , which is estimated as a weighted average across enzyme species.Despite these assumptions and limitations, we believe this work is a step forward in modeling enzyme kinetics.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Figure 1 .
Figure 1.Capturing enzyme heterogeneity Left: several mechanisms that can lead to enzyme heterogeneity.The properties of each enzyme, k cat and K D , are variable and can be described by some unknown probability distribution (with probability density p), shown at the top.Right: flowchart for testing for homogeneity.The end result is a posterior probability of either model (homogeneous and heterogeneous) being the correct model for the dataset, as well as parameter estimates.

Figure 2 .
Figure2.Characterizing the effect of s K on the heterogeneous model Top row: varying s K while keeping all other parameters constant changes the Michaelis-Menten curve.This parameter is the standard deviation of a log-normal distribution.The relative difference between upper and lower quantiles q, for varying standard deviations, are shown in the top right figure.Bottom row: three datasets were simulated under the heterogeneous model, with varying s K .The x axes are concentrations a (log-scale) and the y axes are reaction velocities r, which were simulated with random error e. Bayesian inference was run on these datasets, with both models (homogeneous and heterogeneous).In the second two cases (s K = 2; 3), the heterogeneous model is best, because pðI m = 1Þ is large.However, when s K = 1, the homogeneous and heterogeneous model could not be discriminated between, and therefore the model averaging favored the simpler model I m = 0.

Figure 3 .
Figure3.Well-calibrated simulation study (heterogeneous Michaelis-Menten model) Two hundred datasets were simulated, and the parameters used to generate them were recovered with Bayesian MCMC.The coverage is the percentage of simulations where the true parameter lies in its 95% credible interval (blue lines), compared to the times when its estimate is wrong (black lines).The span is the relative difference between the upper and lower 95% interval points, geometrically-averaged across all replicates and rounded to 2 significant figures.These results show that coverage is close to 95% and that the spans are within an order of magnitude, thus providing confidence in the ability of the method to recover parameters.s K is not part of the homogeneous model, and therefore its coverage is only calculated when pðI m = 1Þ > 0:95.In the bottom middle panel, the estimated model is indicated if it has more than 0.95 posterior support, while '?' denotes uncertainty (i.e., 0:05 < pðI = 1Þ < 0:95).

Figure 4 .
Figure 4. Well-calibrated simulation study (homogeneous Hill model) Each point represents an MCMC chain performed on one of two hundred datasets simulated under known parameters.h is part of the Hill model, and therefore its coverage is only calculated when pðI h s0Þ > 0:95.The true estimated model is indicated if it has more than 0.95 posterior support, while '?' denotes uncertainty (i.e., 0:05 < pðI h s0Þ < 0:95).

Figure 6 .
Figure 6.Testing biological datasets for homogeneity and their use of the Hill coefficient Observations are denoted by open circles, and each colored line represents a sampled fit, under the posterior distribution.Yellow lines indicate a homogeneousmodel-sampled fit, and blue for heterogeneous.Note that the x axes are on a logarithmic scale.The best model was selected using Bayesian model averaging, and the mean R 2 is specified for convenience.

Figure 7 .
Figure 7. Cartoon summaries of the three putatively heterogeneous systems Top: MnP is oxidized by its substrate H 2 O 2 , resulting in a mixture of enzymes with varying catalytic abilities.49Middle: membrane protein CFTR has heightened ATPase activity when open and phosphorylated, however its varying conformations and varying degrees of phosphorylation may perform this reaction at different effective rates.16Bottom: the fish brain mitochondria, and all their various transport, regulatory, and metabolic processes, consume oxygen and produce reactive oxygen species (ROS) in a non-MM manner.22

TABLE
d RESOURCE AVAILABILITY B Lead contact B Resource availability B Data and code availability