Multisite Enzymes as a Mechanism for Bistability in Reaction Networks

: Here, we focus on a common class of enzymes that have multiple substrate binding sites (multisite enzymes) and analyze their capacity to generate bistable dynamics in the reaction networks that they are embedded in. These networks include both substrate − product − substrate cycles and substrate-to-product conversion with subsequent product consumption. Using mathematical techniques, we show that the inherent binding and catalysis reactions arising from multiple substrate − enzyme complexes create a potential for bistable dynamics in such reaction networks. We construct a generic model of an enzyme with n -substrate binding sites and derive an analytical solution for the steady-state concentration of all enzyme − substrate complexes. By studying these expressions, we obtain a mechanistic understanding of bistability, derive parameter combinations that guarantee bistability, and show how changing speci ﬁ c enzyme kinetic parameters and enzyme levels can lead to bistability in reaction networks involving multisite enzymes. Thus, the presented ﬁ ndings provide a biochemical and mathematical basis for predicting and engineering bistability in multisite enzymes.


■ INTRODUCTION
Cellular reaction networks enable cells to remain out of thermodynamic equilibrium and to respond to external cues.The dynamics of these networks enable cellular homeostasis and decision-making. 1,2Many decision-making processes involve so-called bistable dynamics, in which a system can attain two different steady states depending on initial conditions.Bistability is implicated in many cellular decision processes, including cell cycle control, 3 lysis−lysogeny decision, 4 metabolic shifting, 5−7 and persister formation. 8anifestation of bistability requires some mechanism of feedback. 9,10In the case of enzymatic reaction systems, feedback dynamics can arise from transcriptional, substrateor product-based regulation, or via post-translational modification of enzymes.−13 In the case of substrate-and product-based regulation of enzymes, a commonly used model considers an enzyme with two binding sites, where the binding of substrate at one side leads to catalysis, while binding of the substrate or product on the other site alters catalytic rate.In such a two-site enzyme model, both bistability and oscillations are attainable depending on the specific binding mechanisms and the assumed functional forms of the rate equations. 2,11,14,15Despite this wide application of the two-site enzyme model, it is currently not clear how exactly a multisite enzyme facilitates bistability and under which parameter regions and biochemical conditions it does so.This is a relevant question, considering that many enzymes found in central metabolism and signaling pathways are multimers with multiple substrate binding sites. 16Specific examples include dehydrogenases with key metabolic substrates (e.g., phosphoglycerate, malate, and lactate), which are commonly composed of dimers or tetramers with multiple binding sites, 17 and kinases such as phosphofructokinase, which have multiple active binding sites. 18A better understanding of reaction dynamics of multisite enzymes can allow us to predict which naturally existing enzymes might be implementing bistability for cellular decision-making or might be suited for the engineering of bistability through synthetic biology.
Here, we undertake an extensive theoretical study of a generalized model of an enzyme with n-substrate binding sites to derive both a biochemical intuition and a set of mathematical conditions on kinetic parameters for bistability.We use primarily analytical approaches to show that the multisite nature of an enzyme inherently results in a potential for bistability.We then use this insight to derive conditions on the kinetic rate parameters of simple reaction networks with multisite enzymes, which guarantee bistability for some concentration of substrate and enzyme.These findings allow us to predict and outline enzyme engineering strategies that can be employed to achieve bistability in simple reaction networks.

■ RESULTS
To better understand how a multisite enzyme can lead to bistability, we first create a generic model of substrate-(S)-toproduct (P) conversion mediated by an enzyme (E) that has nsubstrate binding sites (Figure 1A).In this initial model, we assume that the total concentration of substrate and product and the total concentration of free and substrate-bound enzymes are conserved (see the Methods section and Supporting Information (SI)).The former assumption is directly applicable when the substrate is a conserved moiety, such as enzyme cofactors or energy and reducing power equivalents (e.g., ATP−ADP and NADH−NAD + pairs). 2,19his assumption is useful to illustrate our results, and relaxing itas discussed belowshows that our main conclusions remain intact for the cases where substrate concentration is freely changing (e.g., through fluxes by other reactions).The latter assumption of the total enzyme concentration being conserved reflects the fact that the time scales of enzyme expression are in most cases slower than the reaction dynamics.
To make the model as generic as possible, we use mass action kinetics with irreversible enzymatic catalysis and consider substrate molecules binding to the enzyme in any order and also irrespective of how many substrates are already bound.As we show in the SI, more restricted assumptions about substrate binding order or affinity do not alter our main conclusions.To exemplify our modeling approach, in the Methods section, we provide the set of reactions arising from the generic model for a two-site enzyme, i.e., n = 2 (see also Figure 1B).For our general n-site model, the full set of reactions can be formally written as where a, b, and c are the kinetic rate constants associated with the individual substrate binding sites i, which are numbered 1 through n.The set [n] is the complete set of binding sites [n] = {1, ..., n}, and ES I and ES J are enzyme complexes in which a given number of substrate molecules are bound, respectively, to a set of sites I and J.In other words, I and J are sets with any number of elements from the list of sites 1 through to n; (I, J ⊆ [n]).For example, for I = {1, 3, 4}, ES I is the enzyme complex where the sites numbered 1, 3, and 4 are bound to substrate molecules (Figure 1A).Additionally, ES J is formed by the binding of a single, additional substrate molecule to ES I , meaning the difference between the sets of I and J in eq 1 is one element.The system defined by eq 1 results in 2 n − 1 enzyme complexes (Figure 1A).Fully Bound and Nonfully Bound Enzyme Complexes Display Distinct Steady-State Dynamics with Increasing Substrate Concentration.We analyzed the above generic model using analytical methods 20 to derive solutions for the steady-state concentrations of all 2 n − 1 enzyme complexes, as functions of the steady-state concentration of substrate ([S]) (see SI for details).We found that the steady-state concentration of any complex (ES I ) is given by Here, E tot is the total enzyme concentration including both free and bound forms of the enzyme, and M = 2 n − 1 − n.The terms |I| and |J| are the number of elements (i.e., bound sites) in a given complex, and thus, the index l, which also appears as an exponent to [S], is over the number of bound substrates.The terms α I,l and α J,l are a function of the kinetic reaction constants associated with each of the enzyme complexes (see SI, Section 1.1 for details).These functions are too complex to be included in open form in eq 2, but for the purposes of this work, it is enough to know that they are polynomials in the kinetic rate constants, with all coefficients positive, thus taking positive values for all parameter values.We note that eq 2 is derived under the most generic case of substrates binding to different enzyme sites in any order; however, we show that eq 2 remains true if we assume more specific binding processes, e.g., binding at a specific enzyme site requiring other sites to be bound with the substrate (see SI, Section 1.1 for details).In such cases, some of α I,l and α J,l might be zero.A close inspection of eq 2 shows that [ES I ] will always be given by a fraction of two polynomials in [S].These polynomials will differ in their degree in [S] unless I is equal to the full set (i.e., I = [n]).This is because, when I ≠ [n], the summations in the denominator and the numerator in eq 2 are over different numbers of bound substrates.Specifically, the summation in the denominator is over all enzyme complexes and the largest degree of this polynomial will be equal to M + n = 2 n − 1, the total number of enzyme complexes.In contrast, the summation in the numerator is over the enzyme complexes that can be generated from the enzyme complex ES I .If ES I is the fully bound enzyme complex, then the degree of the numerator will be equal to that of the denominator, as the largest possible value of the index l would be M If ES I is not the fully bound complex, then the degree of the numerator will be equal to that of the denominator minus the number of empty binding sites in ES I .For instance, if the enzyme has two substrate binding sites, leading to three potential different enzyme complexes (Figure 1B), the degree of the polynomial in the denominator will be three (see SI, Section 1.1).The polynomial in the numerator would have a degree of three for the fully bound complex, while for the two complexes, consisting of one filled and one empty binding sites, it would have a degree of two (see SI, Section 1.1).
The specific structure of eq 2 provides an insight into the behavior of steady-state concentrations of the different enzyme complexes with increasing [S] (Figure 1C).Considering the fact that the polynomials comprising eq 2 have positive coefficients (given by α I,l and α J,l , which are functions of kinetic rates), the steady-state concentration of all enzyme complexes will initially increase from zero with increasing [S].Since eq 2 for the fully bound complex has polynomials of the same degree in the numerator and denominator, the limit value of eq 2 at very high [S] for this complex will be the ratio of the coefficients of the highest degree terms of the numerator and denominator.We show that this ratio is equal to E tot , the total enzyme concentration in the network (see SI, Theorem 1).Thus, for the fully bound enzyme complex, the steady-state concentration will initially increase with increasing [S] and approach finally a positive value given by E tot (Figure 1C, far right).In the case of the nonfully bound enzyme complexes, eq 2 will have a lower-degree polynomial in the numerator than in the denominator, and therefore, its limit value at very high [S] will approach zero.Thus, for the nonfully bound enzyme complexes, their steady-state concentration will initially increase with increasing [S], show at least one peak, and then approach toward zero from above (Figure 1C, far left and middle).
Note that, while Figure 1C shows the behavior of eq 2 for an enzyme with n = 2, the analytical summary presented here is independent of n.It shows that, for a multisite enzyme, we will always have two distinct and qualitatively different curves describing the different enzyme complexes' steady-state concentrations (as exemplified in Figure 1C).From here on, we refer to these two qualitatively distinct types of curves as "positive" and "negative" types.Both positive-and negativetype curves will increase when [S] is small and increasing.At large values of [S], both curves will approach a limit value, with the positive-type curve approaching its limit from below and a negative-type curve approaching its limit from above (Figure 1C).These overall conclusions for curve shapes against small and large values of [S] are independent of the specific values of the kinetic rate parameters.They arise solely because of the polynomial degree structure of eq 2, in other words, from the multisite structure of the enzyme.
The exact shape of the curves for intermediate, increasing values of [S], however, and in particular the number of peaks they will display before approaching the limit value, will depend on the kinetic rate constants of the individual enzyme−substrate complexes (i.e the functions α I,l and α J,l in eq 2).For an enzyme with n = 2, the negative-type curves (of the single-substrate complexes) will always show a single peak and have one inflection point (see SI, Section 1.1).The positive-type curve (of the fully bound, two substrate complex) mostly shows no peaks and is a steady increasing function of [S], but there are kinetic parameters for which it would display peaks, as we discuss below (see SI, Section 2.4).With higher n, the form of the expressions suggests that the negative-and positive-type curves may display multiple peaks, although we have not been able to verify this.Intuitively, and from a biochemical perspective, the positive-type curve can be thought of as a saturation process, in which increasing [S] pushes more enzyme binding sites to be filled, ultimately leading to an increase of the steady-state concentration of the fully bound enzyme complex.Correspondingly, the steady-state concentrations of the nonfully bound enzyme complexes decrease with increasing [S], giving rise to the negative-type curve.
Negative-Type Curves of Nonfully Bound Enzyme Complexes Underpin the Potential for Bistability.We now consider the catalytic flux through each enzyme complex.We refer to the catalytic flux through complex ES I as which is simply the concentration of the complex, [ES I ], multiplied by the sum of the catalytic rate constants of product formation from ES I .Note that the steady-state value of V ES I S→P will be a function of the steady-state complex concentration.Furthermore, the total catalytic flux through the enzyme, V S→P , will be given by the sum of the individual fluxes through each of its complexes.By the analysis above, V S→P tends to E tot , times the sum of the catalytic rate constants of the fully bound complex.To illustrate the idea for general n, we consider first the example case for an enzyme with n = 2 (shown in Figure 1B,C).It is easier to graphically understand how bistability arises in this system if we analyze the behavior of the catalytic fluxes against S sum = S tot − [P], where S tot is a constant describing the combined amount of product and free and bound substrates (see SI, Section 1.1).As we show in the SI, for n = 2, S sum is an increasing function of [S] and hence the qualitative behavior of V S→P against increasing [S] or S sum is the same.
In Figure 2, we show V ES I S→P against S sum for two different parameter sets, and as expected, we see that the behavior of as given by eq 2 and shown in Figure 1.In the example shown in Figure 2A, where we have the same parameters as in Figure 1C, the total catalytic flux V S→P is dominated by the fluxes through the nonfully bound complexes, and as such, V S→P displays a negative-type behavior in S sum .In Figure 2B, we see the results for a second set of parameters, where V S→P is dominated by the flux through the fully bound complex, and as a result, it displays a positive-type curve in S sum .As illustrated by these examples, which type of behavior V S→P displays depends on the specific values of the kinetic rate constants of the individual enzyme complexes.We now consider the shape of the V S→P curve in the context of a reaction system.To start with, we consider a simple scenario, involving a backreaction from the product to substrate, creating a reaction cycle (see Figure 2C).We initially assume that the product-to-substrate conversion is a nonenzymatic, hydrolysis-type reaction, governed by a constant k h (note that, below and in the SI, we relax this assumption without the loss of the presented conclusions).The catalytic flux of this backreaction, V P→S , is given by k h × [P] and therefore behaves linearly with increasing S sum .This linear relation has slope −k h and intercept S tot (Figure 2C).When we plot V S→P and V P→S against S sum on the same plot, the intersection points represent the steady states of the reaction system, i.e., points where the product formation flux, V S→P , equals that of product loss, V P→S .Using the fact that V P→S is a line with a negative slope, we can see that a negativetype V S→P curve opens the possibility to have three intersections between V S→P and V P→S and therefore three steady states.Three steady states are the hallmark of bistability, and indeed, for this parameter set, our model displays bistability, where different starting conditions can lead to different steady-state dynamics (Figure 2D).Since adjusting the value of S tot results in shifting the V P→S line along the xaxis, we can graphically see that as long as k h is below a certain threshold value, there will be some value of S tot that ensures three intersections.In other words, tuning the S tot value would allow shifting the V P→S line across the x-axis on Figure 2C, until three intersections with the V S→P curve are obtained.
While we analyze a system with n = 2 and a sample parameter set in Figure 2, we can use the above discussion to draw a general conclusion that will be true for any n.If V S→P is of the negative type and its slope at the inflection point is smaller than that of V P→S (that is, −k h ), then the curves will intersect three times if the line V P→S passes through the inflection point.The slope of V S→P at the inflection point depends on E tot and the reaction rate constants, while the slope of V P→S at the inflection point depends on k h .This graphical analysis, therefore, provides an intuition about why having a V S→P of the negative type and with a slope at its inflection point smaller than −k h provides a route to bistability in a system with a multisite enzyme for some value of S tot .On the The first model, on the far left, considers a cyclic reaction from substrate to product and back to the product, with the latter reaction being nonenzymatic.The second model, on the middle panel, has the same structure but considers the reverse reaction from product to substrate being enzymatic.The last model, on the far right, considers in-and out-fluxes for the substrate and product, with the cyclic product-to-substrate reaction being optional (including it does not alter conclusions; see SI, Section 2.3).For simplicity, each system is shown with a two-site enzyme model and with only one of the catalytic reactions (on one example binding site).Further model details are provided in the SI.Note that the mathematical analysis presented in the main text considers an n-site model, with all possible binding and catalysis reactions.(B) The core inequality, as shown in eq 3 and common to all of the cases considered, is written for the generic, n-site model.This inequality characterizes when V S→P is of the negative type.We note that the right side of this equation corresponds to only the sum of catalytic rates from the fully bound enzyme complex, as depicted in the cartoon below.The left side of the inequality involves both catalytic rates and equilibrium constants of those enzyme complexes that are unbound only on one site.contrary, if V S→P is of the positive type and does not display any peaks (as shown in Figure 2B), bistability is precluded as V P→S cannot intersect V S→P in more than one point.When V S→P is of the positive type and displays a single or multiple peak, there is again the possibility for three intersection points and bistability (see SI, Section 2.4).In summary, this graphical discussion shows that a negative-type curve for V S→P guarantees three steady states after appropriately choosing the other relevant parameters (e.g., k h and S tot ).
Kinetic Rate Parameter Conditions That Guarantee Multiple Steady States in a Reaction System with a Multisite Enzyme.To formalize and generalize the graphical considerations made above, we take a mathematical approach to determining conditions on kinetic parameters that result in multiple steady states.The idea is to identify the conditions when V S→P is of the negative type, that is, when it converges to its limiting value from above, and use these conditions to guarantee that V P→S and V S→P will intersect at multiple points.
We find that V S→P is of the negative type exactly when the following condition holds (SI, Section 1.2) Here, = represents the catalytic rate constants of the fully bound enzyme complex, where catalysis occurs at the i-th binding site (see also Figure 3).We note that the condition defined by eq 3 is aligned with the graphical analyses discussed in the previous sections (Figures 1 and 2).There, we have shown that the curve type of V S→P is determined by whether the fully bound or nonfully bound enzyme complexes are dominating the dynamics of catalysis.In line with these arguments, for eq 3 to hold and hence for V S→P to be of the negative type, the sum of the catalytic rate constants for the n − 1 nonfully bound complexes, each adjusted by the contribution of that complex in the system dynamics (represented by their K m 's), has to be greater than that of the catalytic rate constants of the fully bound complex.
To reiterate, eq 3 determines the condition for V S→P to be of the negative type.For bistability to arise, the negative-type V S→P needs to be complemented with an appropriate V P→S curve (or line) and this latter aspect will depend on the system in which the multisite enzyme is embedded.We first study the simple case of a nonenzymatic backreaction from product to substrate (see the next section for results of alternative reaction systems).We find that we are guaranteed to have three positive steady states in such a cyclic reaction system, for some values of S tot and E tot , if the reaction rate constants satisfy the following condition (SI, Section 2.1) This condition is identical to eq 3 but with the extra term k h on the right-hand side.This is again in line with our analysis above.First, if eq 4 holds, then eq 3 must also hold, and hence, V S→P is of the negative type.Second, eq 4 tells us that k h cannot be larger than a certain amount.This ensures that, with the appropriate choice of S tot and E tot , V P→S passes through the last inflection point of V S→P with a slope larger than the slope of V S→P at that point.As discussed, this gives rise to multiple steady states.
Conditions for Multiple Steady States Exist for Different Reaction Systems Involving a Multisite Enzyme.Using the same approach as above, we expanded our analysis to other realistic reaction motifs featuring a multisite enzyme.We considered two common motifs, involving a nonenzymatic or enzymatic backreaction from the product to substrate (Figure 3A, far-left and middle panels), or in-and out-fluxes of both substrate and product (Figure 3A, far-right panel).The former case represents two enzymes creating a cyclic reaction motif and is commonly found in metabolism and in signaling systems. 2,14,15,19The latter case represents another widely applicable scenario, where any upstream and downstream reactions can generate or consume the substrate and product.In this case, there is no assumption of the total substrate amount being conserved.
For each of the cases depicted in Figure 3A, we found that the existence of multiple steady states is guaranteed by an inequality almost identical to eq 4 (see SI Sections 2.2 and 2.3).In the case of a reaction system with an enzymatic backreaction from the product to substrate, the catalytic rate constant of the backreaction replaces k h in eq 4. In the case of the reaction system involving fluxes of the substrate and product, k h is eliminated entirely from the inequality, that is, the inequality reduces simply to eq 3.These resulting inequalities need to be supplemented with a distinct choice of additional parameters.For the system with enzymatic backreaction, eq 4 guarantees multiple steady states after appropriately selecting S tot , E tot , and the conserved total amount of the enzyme catalyzing the backreaction from the product to substrate (SI, Section 2.2).For the system with fluxes, eq 3 guarantees multiple steady states after appropriately choosing E tot and flux rate constants (SI, Section 2.3).So, in this case, the possibility of multiple steady states is not conditioned on the value of S tot as the total amount of substrate is no longer conserved.
The key intuitive message, as depicted in Figure 3B, is that a key sufficient mechanism for the existence of multiple steady states is related to the dynamics of two distinct sets of enzyme complexes, those that are fully bound and those that have one binding site empty.When the kinetics of the latter dominates over that of the former and eq 3 is satisfied, a negative-type V S→P curve emerges from the multisite enzyme dynamics and multiple steady states are guaranteed to exist in some parameter regime in the system.
It is important to note that especially with increasing n many multiple steady states may arise and not just three.We note that a formal analysis of the stability of each steady state cannot be done using the presented general framework.In the case of systems with n = 2 and 3, we have sampled kinetic parameter values satisfying eq 4 and found that when the system displays three steady states, then bistability arises, showing that at least two steady states are stable.Finally, we also note that eqs 3 and 4 do not define necessary conditions for multiple steady states but rather conditions that guarantee multiple steady states.As we argued above, there can be parameter sets that lead to a positive-type V S→P curve with multiple peaks and therefore still lead to multiple steady states without fulfilling eq 4 (SI, Section 2.4).
Enzyme Parameters in the Physiological Ranges That Satisfy eq 3 Permit Bistability.As described above, eqs 3 and eq 4 describe conditions on the catalytic rates and K m constants that are guaranteed to result in multiple steady states for some set of S tot and E tot values.To identify ranges of these latter parameters, we used numerical and analytical methods with the 2-site enzyme model with a cyclic reaction motif involving a nonenzymatic backreaction as a case study (first panel of Figure 3A).We have chosen kinetic parameters in a physiological range using available information from the literature on multisite enzymes involved in cyclic reaction systems (see Methods).We then derived a bifurcation diagram for the parameters S tot and E tot (see Methods).We find that for physiologically relevant kinetic parameters, there is a relatively wide range of S tot and E tot values allowing for multiple steady states, but E tot is always much smaller than S tot (Figure 4A, red area bounded by dashed lines).In other words, the manifestation of multiple steady states in this cyclic reaction scheme happens in a regime of substrate-saturated enzymes.In fact, for this reaction system, we find that the relation S tot > n • E tot needs to hold for systems satisfying eq 4 to display multiple steady states (see SI, Section 2.1).
How would changing kinetic parameters affect the S tot and E tot ranges permitting multiple steady states?As discussed above, S tot determines the intersection point of the V P→S line with the x-axis, while E tot determines the height of the V S→P curve.We can therefore expect that kinetic parameters affecting the slope and shape of the V P→S line and the V S→P curve will alter the S tot and E tot ranges permitting multiple steady states.In line with this prediction, we find that decreasing k h and increasing the catalytic rates of the nonfully bound enzyme complexes widen the S tot and E tot ranges for multiple steady states (Figure 4A; regions bounded by straight and dotted lines).The latter creates this effect by changing the slope of the V P→S line, while the latter by changing the height of the V S→P curve.
In the case of the reaction system with substrate and product fluxes (Figure 3A, left-most panel), i.e., where S tot is not constant, the bistable regime is determined by enzyme kinetic parameters, substrate in-and out-flux, product out-flux, and E tot (see SI Section 2.3).For this case, we derived a bifurcation diagram for substrate in-flux and product out-flux rates for a given physiologically realistic E tot and found that changing E tot can result in widening of the bistable regime for these two parameters (Figure 4B).
Engineering Routes toward Implementing Bistability in Multisite Enzymes.We have provided in eq 4 (see also eq For the region bounded by the straight lines (covering all of the gray area and some of the red area), the only parameter altered was the hydrolysis rate of the product; k h =10 2 min −1 .For the region bounded by the dotted lines (covering all of the blue and gray areas, and some of the red area), the two parameters altered were the hydrolysis rate of the product and the catalytic rate of one of the single-bound complex; k h = 10 2 min −1 and c 0,1 = 10 6 .Note that the left boundary of the regions bounded by the straight and dotted lines overlap.(B) Two-parameter bifurcation diagram for the reaction system with free substrate and product fluxes (as shown on the right-most cartoon on Figure 3A) and involving a 2-site enzyme.The diagram shows the regime with three steady states for varying substrate in-flux rate k S,in (y-axis) and product out-flux rate, k P,out (x-axis).Parameters used were as for the straight-line case of part A, and with additional parameters set as; k S,out = 10 min −1 , k P,in = 0 (no product in-flux).The parameter E tot was set to 4.15 × 10 −5 and 10 −4 M for the areas bounded by the straight and dashed lines respectively.
7 and related equations in the SI) mathematical criteria for a multisite enzyme to display multiple steady states when complemented by appropriate choices of total enzyme and substrate concentrations.One of the key findings from these analyses is that having a lower catalytic rate from the fully bound enzyme complex (the term c [n]\ {i},[n] ) in eqs 3 and 4, compared to that of nonfully bound complexes, favors bistability.Since substrate inhibition could be seen as a manifestation of the fully bound enzyme complex displaying lower catalytic activity, enzymes displaying substrate inhibition could present possible candidates for testing the presented theory.In Table 1, we list a set of multisite enzymes that are shown to display substrate inhibition and note that these enzymes might already be satisfying the conditions presented in eqs 3 and 4. Indeed, for two of the listed enzymes, namely, phosphofructokinase and lactate dehydrogenase, in vitro reconstitution experiments implementing cyclic reaction systems and using experimental conditions of substratesaturated enzymes have directly demonstrated the existence of bistability. 21,22As utilized in those studies, experimental in vitro studies can directly look for bistability in all of the listed, and similar, multisite enzymes.
Systems not displaying bistability can be further manipulated toward changing the K and c terms as presented in eqs 3 and 4 (or eqs 7 and 8) to fulfill the mathematical criteria given in these equations.To this end, mutations altering either the Michaelis−Menten or catalytic constants (i.e., the K and c terms in eqs 3 and 4 or eqs 7 and 8) have been described previously for several multisite enzymes. 18,23As an alternative to mutations, one could also use inhibitory compounds to specifically affect the catalytic rate of the fully bound enzyme complex (term c [n]\ {i},[n] in eqs 3 and 4).In this context, inhibitors have been identified that specifically alter catalytic activity, rather than the Michaelis−Menten constant, for some multisite enzymes such as 3-phosphoglycerate dehydrogenase. 18In that case, the inhibition mechanism is related to the slowing (or complete stopping) of the movement of substrate binding protein domains required for catalysis, 18 but it remains to be tested if such an inhibition mechanism can be made to specifically alter the catalytic rate of the fully bound enzyme complex.As a third possible route for engineering, domainbased protein engineering (e.g., as in ref 24) could potentially be used to add further substrate binding sites without catalytic activity on a single-substrate binding enzyme.
Another experimental direction for engineering bistability would be to focus on altering enzyme and substrate concentrations for a given system featuring a multisite enzyme that is already shown to display low activity at high substrate binding, i.e., substrate inhibition 25 (see also Table 1).As we show in Figure 3, even for enzymes satisfying conditions given in eq 4, their bistability will be manifest only for certain ratios of substrate-to-enzyme concentration and particularly when S tot > n • E tot is satisfied, as discussed above.Thus, experimentally manipulating S tot and E tot might provide means to implement bistability in an in vivo or constituted enzymatic system.To this end, the use of quantitative control of gene expression levels based on CRISPR 26 can be explored to control enzyme levels directly.
A multisite enzyme with bistable dynamics will display bimodal enzymatic fluxes, which can be seen through the measurement of substrate and product levels under different starting conditions (as shown in Figure 2).This approach is used to demonstrate bistability in lactate dehydrogenase and phosphofructokinase enzymes, which are reconstituted in vitro and their cosubstrate, substrate, and product levels are monitored spectroscopically or chemically. 21,22We propose combining such an in vitro assay with the engineering strategies described above to specifically engineer bistability in a multisite enzyme.An ideal case study for such engineering could be the isocitrate dehydrogenase, which is already identified with potential for bistability according to the presented theory due to its observed substrate inhibition and multisite nature (see Table 1) and is postulated to form a bistable reaction system. 15his enzyme (and similar enzymes as listed in Table 1) can be incorporated into a reconstituted reaction system in vitro, and the impact of different mutations (generated through random or targeted mutagenesis), inhibitory compounds, and enzyme concentration on their bistability can be directly evaluated by monitoring enzymatic fluxes under different initial conditions.This approach should allow the identification of bistability and bistability-causing mutations.The linkage to the presented theory can be further made by measuring kinetic parameters of mutant and wild-type enzymes resulting from such studies.Any such engineered bistable enzymes, or naturally existing ones, can be utilized to study the effect of bistability on cell physiology and for implementing dynamical control of pathway fluxes.The latter route is important for addressing the identified key challenge of dynamic control of metabolically engineered pathways, 27,28 as well as recent aims of decoupling cellular production and growth stages through bistable systems. 29In the case of isocitrate dehydrogenase, which is shown to be involved in the metabolic engineering of citrate overproduction, 30 the engineering of bistability can allow cells to switch between a low-and high-flux regime for isocitrate dehydrogenase flux depending on growth conditions (e.g., glucose levels), thereby providing the possible control of downstream citrate production.More broadly, natural or mutant enzymes with demonstrated bistability can be utilized to study physiological consequences of bistability such as the generation of metabolically distinct subpopulations. 31METHODS Core Biochemical Model.We considered first a core model involving an enzyme with multiple substrate binding sites, each able to convert the substrate into a product, as shown in Figure 1.For this model, we assumed that the total enzyme concentration and the total substrate concentration, which is the free substrate, the substrate bound to enzyme, and the product, are conserved.We relaxed the latter assumption in subsequent models that were built from this core model.For the core model, the resulting binding and catalytic reactions for an enzyme with n-binding sites are given in eq 1.Additional reactions in the subsequent models and involving the product, and sometimes the substrate, are considered, either as occurring with a constant rate or mediated by an additional enzyme.Our mathematical analyses consisted of writing ordinary differential equations (ODEs) for such reaction systems using mass action kinetics.The ODEs for the core, general model shown in Figure 1, as well as the alternative models shown in Figure 3, are provided in full in the SI along with the detailed derivations leading to eqs 2−4.As an illustration, we provide here the reaction system for the core model, for n = 2, i.e., a two-site enzyme where the single-and double-bound enzyme complexes are denoted as ES 1 , ES 2 , and ES 1,2 , respectively.The corresponding set of ODEs resulting from this reaction system can be written using mass action kinetics for each of the reactions shown in eq 5, as we have done in the provided MATLAB code (see SI file1).The conservation relations for this system are Note that there are three enzyme complexes for this system; thus, 2 n − 1 = 3.The specific form of eq 4 for this system is thus given by where K 1 and K 2 are the Michaelis−Menten constants for enzyme complexes ES 1 and ES 2 featured in eq 5. Writing these Michaelis−Menten constants in their open form (using the kinetic rate constants shown in eq 5) and after algebraic manipulation, eq 7 is equivalent to As discussed in the main text, these inequalities show the key requirement on specific kinetic parameters for bistability to occur.We can see from eq 7 that the requirement for bistability is for the sum of catalytic rates from single-bound enzyme complexes, adjusted by their Michaelis−Menten constants, to be larger than the sum of the catalytic rates from fully bound enzyme complexes plus the kinetic constant for hydrolysis.
Symbolic and Numerical Computations.For all symbolic computations, utilized in finding steady-state solutions and deriving mathematical conditions on rate parameters, we used the software Maple 2020.For simulations, run to numerically analyze selected systems, we again used Maple, or the MATLAB package, with the standard solver functions.
Bifurcation Analysis and Physiologically Realistic Kinetic Parameters and S tot and E tot Ranges.To analyze if multiple steady states would be realized in physiologically realistic parameter regimes, we used a cyclic reaction system with a two-binding-site enzyme (Figure 4A).For such an enzyme, we have used kinetic parameter values in physiologically feasible ranges as found in the literature and listed below. 16,32,33We then used our mathematical condition shown in eq 4 and bifurcation analyses to derive the S tot and E tot ranges that guarantee multiple steady states.The analysis was performed using cylindrical algebraic decomposition in Maple, using the package RootFinding[parametric]. 34 As an example, the kinetic rate values used for Figure 2, as listed on its legend, satisfy eq 4 and hence would result in multiple steady states when combined with any combination of S tot and E tot that are in the permissible range shown in Figure 4. Physiologically realistic kinetic parameter ranges (based on 16,32,33) considered were 10 7 −10 10 M −1 min −1 for substrate−enzyme binding, 10 2 −10 6 min −1 for substrate dissociation from a substrate−enzyme complex, 50−10 7 min −1 for catalytic rates of enzyme complexes and hydrolysis rates (i.e., k h ), 10 −6 −10 −2 M for K M values, and 10 −6 −10 −2 M and 10 −8 −10 −4 M for S tot and E tot , respectively.

■ DISCUSSION
We have shown that multisubstrate binding enzymes have an inherent capacity to generate bistability when placed within a reaction system.Specifically, the very act of an enzyme binding two or more molecules of the same substrate is guaranteed to result in a specific nonlinear relation between the substrate amount and catalytic flux rate (V S→P ) in a certain parameter regime (we called the resulting relation a negative-type curve in the main text).When the multisubstrate−enzyme is placed within a reaction system, this inherent dynamical feature of a negative-type curve then guarantees the emergence of multiple steady states.The wider reaction systems embedding a multisite enzyme can involve either substrate−product− substrate cycles or systems involving open substrate and product fluxes.
These types of reaction systems, as well as multisite enzymes embedded in them, are common occurrences in metabolic and signaling pathways.Dehydrogenases and kinases, for example, are commonly involved in substrate-to-product cycles (as shown in Figure 3A, middle panel), via either redox cycling or phosphorylation/dephosphorylation of substrate−product pairs.Examples include reactions involving dehydrogenases such as lactate or glutamate dehydrogenase, 21 and kinase/ phosphatase pairs such as those involved in the conversion of fructose-6-phosphate. 22,35,36 The case with substrate and product fluxes (Figure 3A, far-right panel) is a particularly generic scenario, where there is no mass conservation assumption with regard to the substrate and product and no requirement for a cyclic reaction motif.In these different, common reaction systems, we demonstrate that a multisite enzyme can lead to bistable dynamics.This is because the negative-type V S→P curve is an inherent feature of the multisite enzyme and therefore independent of downstream product (and substrate) conversions.Thus, any arrangement of a reaction system resulting in substrate and product conversion dynamics that is capable of intersecting a V S→P curve of a negative type three times will result in a system capable of multiple steady states, as we show here.
To directly ascertain bistable parameter regimes, we derived here a mathematical inequality (eq 3) that guarantees the V S→P to be of the negative type.This inequality constitutes the core part of additional inequalities (see eq 4 and SI) that are derived for different, and common, scenarios of reaction systems embedding a multisite enzyme and that guarantee the existence of multiple steady states in them.Two key, biochemical intuitions arise from these mathematical inequalities.The first is that the different enzyme complexes arising in a multisite enzyme, which have a different number of substrates bound provide a key role for the potential for bistability.These complexes, even if they have fast substrate binding dynamics, provide a form of hidden feedback in the reaction system, which is an essential feature for bistability to arise.In this context, we note that cooperative binding, i.e., higher binding rates for subsequent substrates, is not part of the mathematical requirements given in eqs 2 and 4. Thus, the presented theory shows that multisite enzymes can display bistability even without cooperative binding dynamics.
The second key biochemical intuition arising from the presented theory is that bistability within a system containing a multisite enzyme requires nonfully bound enzyme complexes to "outcompete" the fully bound complex in terms of catalysis (or flux) from the substrate to product.This relates our work to the concept of "substrate inhibition", which is observed in the case of many multisite enzymes and specifically dehydrogenases and kinases 25 and which is commonly attributed to allosteric effects (i.e., substrate binding also at a noncatalytic, regulatory site on the enzyme).In our case, we emphasize that we do not consider extreme allosteric effects such as binding at one site completely abolishing catalytic activity; however, we note that the general model we describe here would have in it nested all types of different allosteric models (e.g., those assuming no catalytic activity for some enzyme−substrate complexes or those assuming a reduction in catalytic rate with increasing substrate concentration, i.e., substrate inhibition).Indeed, when the criterion on kinetic parameters given in eq 3 is satisfied, the resulting dynamics of catalysis rate with increasing substrate concentration (as shown in Figure 2A) will be similar as seen with allosteric models or models implementing substrate inhibition.Whether, in the case of specific, natural enzymes displaying substrate inhibition, the fully bound enzyme complexes have indeed specifically lower catalytic rates than complexes with nonfully bound complexes, needs to be determined through enzyme kinetics experiments.
In addition to the presented inequalities being satisfied, bistability also requires additional system parameters to be chosen appropriately when we consider systems embedding multisite enzymes.We find that these additional system parameters, total substrate and enzyme concentrations, as well as kinetic rate constants of additional reactions leading to multiple steady states, exist within physiologically feasible parameter values obtained from enzymatic studies.A key aspect that we note, in the case of the cyclic system, is that the total substrate levels (i.e., substrate and product combined) need to be larger than the total enzyme concentration.This condition is found to be satisfied for many enzymes in vivo. 32,33n line with these findings demonstrating physiological feasibility, bistability in systems involving cyclic reaction motifs is observed when multisite enzymes are reconstituted in vitro; for example, using pyruvate kinase, lactate dehydrogenase, or isocitrate dehydrogenase enzymes and their corresponding partners, bistability has been demonstrated experimentally. 15,21,22In the case of systems with open substrate and product fluxes, eq 3 guarantees multiple steady states after appropriately choosing E tot and flux rate constants.Interestingly, in this case we find that the tuning of total enzyme levels, which can be implemented with gene expression control, can widen or limit the bistable parameter regime.Therefore, our findings of bistability and the parameter regimes it is manifested in can be of wide relevance for the study of a large range of cellular reaction systems.
Reaction system dynamics are implicated to possess a level of autonomous regulation. 1,2Our findings show that multisite enzymes can indeed provide such regulation by providing reaction systems with the capability of bistability.When bistability is realized, this will manifest itself as two different steady-state concentrations, among which the system can switch.Multisite enzymes can provide a simple mechanism to achieve such higher-level functions.As such, they are complementary and can work in conjunction with other mechanisms that are found to lead to nonlinear dynamics in enzyme reaction systems.−40 Enzyme sequestration involves the assumption of noncatalytic enzyme−substrate complexes, and as such, it is nested in the generic multisite model presented here.We particularly note that when an enzyme sequestration model assumes the noncatalytic complexes to be the fully bound enzyme complexes, then our theory shows that multiple steady states will be guaranteed.Multisite phosphorylation, as studied in the literature, is assumed to involve single-site enzymes acting on substrates with multiple post-translational sites.Thus, it is conceptually similar to the presented theory, in the sense that it highlights the role of multiple enzyme−substrate complexes in generating multistability.We note that it is possible that the two mechanisms can be combined in natural or engineered systems.In other words, we can imagine multisite enzymes catalyzing the conversion of substrates that have multiple modification sites.Such a combination of both scenarios should increase the possibility of the system to display multiple steady states or different dynamical behaviors such as oscillations.
The presented findings provide clear experimental routes toward generating or removing bistability in natural reaction systems or engineered enzymes through the control of kinetic parameters or expression levels with synthetic biological approaches. 24,27The engineering principles described here for bistability can be further extended to explore the possible sources of multistability and oscillatory dynamics, both of which are observed in models with multisite enzymes with flux, 2,14,41

Figure 1 .
Figure 1.(A) Cartoon representation of a generic n-site model, where n-E indicates an enzyme with n-substrate binding sites.The substrate binding sites are numbered in a consecutive fashion, and substrate-bound sites are shown in blue.Note that there are 2n − 1 possible substrate− enzyme complexes.(B) Cartoon representation of a two-site enzyme model.The substrate (S) and product (P) are shown in blue and red, respectively.Substrate binding is allowed in any order on each site, and both sites are assumed to have catalytic activity.The three possible substrate−enzyme complexes are shown on the right (see the Methods section for reactions and differential equations for this two-site enzyme model).(C) The steady-state concentration of each of the substrate−enzyme complexes with increasing concentration of substrate.The parameters, as listed in eq 5, are set to the following values for these simulations: a 0,1 = a 0,2 = a 1,{1,2} = a 2,{1,2} = 10 8 M −1 min −1 , b 0,1 = b 0,2 = b 1,{1,2} = b 2,{1,2} = 10 4 min −1 , c 0,1 = 10 5 min −1 , c 0,2 = 1.5 × 10 5 min −1 , c 1,{1,2} = c 2,{1,2} = 10 3 min −1 , S tot = 2.31 × 10 −3 M, and E tot = 4.15 × 10 −5 M. Panels from left to right show the steady-state concentrations of the two single-substrate complexes and the fully bound complex.On the right-most panel, the dashed line indicates total enzyme concentration.

Figure 3 .
Figure 3. (A) Three different reaction systems, embedding a multisite enzyme, considered in this work.The first model, on the far left, considers a cyclic reaction from substrate to product and back to the product, with the latter reaction being nonenzymatic.The second model, on the middle panel, has the same structure but considers the reverse reaction from product to substrate being enzymatic.The last model, on the far right, considers in-and out-fluxes for the substrate and product, with the cyclic product-to-substrate reaction being optional (including it does not alter conclusions; see SI, Section 2.3).For simplicity, each system is shown with a two-site enzyme model and with only one of the catalytic reactions (on one example binding site).Further model details are provided in the SI.Note that the mathematical analysis presented in the main text considers an n-site model, with all possible binding and catalysis reactions.(B) The core inequality, as shown in eq 3 and common to all of the cases considered, is written for the generic, n-site model.This inequality characterizes when V S→P is of the negative type.We note that the right side of this equation corresponds to only the sum of catalytic rates from the fully bound enzyme complex, as depicted in the cartoon below.The left side of the inequality involves both catalytic rates and equilibrium constants of those enzyme complexes that are unbound only on one site.
−Menten constant of the enzymatic conversion with enzyme complex ES I (as in eq 1), for the enzyme complexes with n − 1 sites bound, and c I\{i},I are the catalytic rate constants as in eq 1.The term c [n]\ {i},[n]

Figure 4 .
Figure 4. (A) Two-parameter bifurcation diagram for the reaction system shown on the inset of Figure 2C and involving a 2-site enzyme with a P to S backreaction.The diagram shows the regime with three steady states for varying S tot (x-axis) and E tot (y-axis) values (in M) for three different sets of physiologically relevant enzyme kinetic parameters.The kinetic parameters used for the region bounded by the dashed lines (covering all of the red area) were:a 0,1 = a 0,2 = a 1,{1,2} = a 2,{1,2} = 10 8 M −1 min −1 , b 0,1 = b 0,2 = b 1,{1,2} = b 2,{1,2} = 10 4 min −1 , c 0,1 = 10 5 min −1 , c 0,2 = 1.5 × 10 5 min −1 , c 1,{1,2} = c 2,{1,2} = 10 3 min −1 , k h = 0.5•10 3 min −1 .For the region bounded by the straight lines (covering all of the gray area and some of the red area), the only parameter altered was the hydrolysis rate of the product; k h =10 2 min −1 .For the region bounded by the dotted lines (covering all of the blue and gray areas, and some of the red area), the two parameters altered were the hydrolysis rate of the product and the catalytic rate of one of the single-bound complex; k h = 10 2 min −1 and c 0,1 = 10 6 .Note that the left boundary of the regions bounded by the straight and dotted lines overlap.(B) Two-parameter bifurcation diagram for the reaction system with free substrate and product fluxes (as shown on the right-most cartoon on Figure3A) and involving a 2-site enzyme.The diagram shows the regime with three steady states for varying substrate in-flux rate k S,in (y-axis) and product out-flux rate, k P,out (x-axis).Parameters used were as for the straight-line case of part A, and with additional parameters set as; k S,out = 10 min −1 , k P,in = 0 (no product in-flux).The parameter E tot was set to 4.15 × 10 −5 and 10 −4 M for the areas bounded by the straight and dashed lines respectively.

Table 1 .
Selected Multisite Enzymes That Are Experimentally Shown to Exhibit Substrate Inhibition a a Only substrates displaying substrate inhibition are listed.