Effect of magnitude and variability of energy of activation in multisite ultrasensitive biochemical processes

Protein activity is often regulated by ligand binding or by post-translational modifications such as phosphorylation. Moreover, proteins that are regulated in this way often contain multiple ligand binding sites or modification sites, which can operate to create an ultrasensitive dose response. Here, we consider the contribution of the individual modification/binding sites to the activation process, and how their individual values affect the ultrasensitive behavior of the overall system. We use a generalized Monod-Wyman-Changeux (MWC) model that allows for variable conformational free energy contributions from distinct sites, and associate a so-called activation parameter to each site. Our analysis shows that the ultrasensitivity generally increases as the conformational free energy contribution from one or more sites is strengthened. Furthermore, ultrasensitivity depends on the mean of the activation parameters and not on their variability. In some cases, we find that the best way to maximize ultrasensitivity is to make the contribution from all sites as strong as possible. These results provide insights into the performance objectives of multiple modification/binding sites and thus help gain a greater understanding of signaling and its role in diseases.

Multisite protein modification is ubiquitous in gene regulation and signal transduction, often in the form of multisite phosphorylation. Many models of multisite ultrasensitivity are available in the literature, but they usually assume that all sites contribute equally to the activation of the multisite target. In this work, we relax this assumption and carry out computational and mathematical analysis of a multisite system in which the conformational free energy contribution varies across sites. We find that the ultrasensitivity of the system tends to increase (with some exceptions) when the conformational free energy contributed by any given site is strengthened. Our analysis predicts that all active sites should have approximately the same conformational free energy contribution, a property observed in proteins with unstructured modification domains and bulk electrostatics. We were also able to predict from first principles an energy range of -2 to -4 kcal/mol per site that effectively maximizes ultrasensitive behavior. This prediction is consistent with Introduction Cellular systems rely heavily on signal transduction and environmental sensing pathways to successfully respond to internal and external environmental signals and conditions. To regulate signal transduction cascades, mammalian cells use ligand binding or post-translational modifications (PTMs) such as protein phosphorylation, methylation or ubiquitination. Several forms of disease can arise when there are defects in signal transduction pathways, including cancer, diabetes, and heart disease. Many proteins regulated by ligands or PTMs are multisite proteins, that is, they have multiple sites on which they can be modified or where a ligand can bind. For example, activation of mitogen activated protein kinases requires phosphorylation on two sites [1], and the hemoglobin tetramer has four sites where oxygen can bind ( [2] and references therein). In fact, some proteins have more than 150 modification sites [3]. Ligand binding/PTMs can either promote or inhibit protein activity through conformational changes [4,5], and can influence the target's enzymatic activity, location, stability, or interactions with other macromolecules [6].
A common role for multisite modifications lies in the creation of switch-like, or ultrasensitive, dose response curves [7][8][9]. These are positive, monotonically increasing, sigmoidal functions that have two important properties: first, they respond minimally to low levels of input; second, once the input is sufficiently large, they switch from a low output to a near maximal output in response to a relatively small increase in the input. In other words, ultrasensitive systems can both filter out low-level noise and respond with a high gain over an appropriate range of input. Ultrasensitivity has important roles in signal transduction, and a widely-studied problem is how to implement such dose responses using common biochemical reactions [2]. A classical 1965 model by Monod, Wyman and Changeux (MWC) [10] uses multisite modifications to create ultrasensitive responses. This model remains highly influential today [11][12][13][14]. In the MWC model, the target molecule/receptor can be in either an active (relaxed) conformation or an inactive (tense) conformation, and ligand binding/PTMs can influence the probability that the target is in one state or the other. One way to envision this is that ligand binding promotes a conformational change that flips the target from inactive to active (or from active to inactive in the case of an inhibitory ligand). In the MWC model, this is equivalent to the point of view that the ligand binds preferentially to the active conformation, a phenomenon known as conformational selection. In multisite MWC models, there are multiple binding sites for ligand, each of which can be either empty or bound. Such models exhibit cooperativity in ligand binding, as the binding of some ligands to the target will promote flipping to the active state, and in the active state, all binding sites have a higher affinity for ligand. In other words, the presence of ligand increases the probability of the receptor existing in the state with higher ligand affinity, thereby increasing the probability of the next ligand binding. In addition to 'cooperativity', the term 'allostery' is frequently used in conjunction with MWC models, and refers to the effect that one ligand binding to the target has on additional "distant" binding sites in the same molecule, as well as to the effect that ligand binding has on the conformational change that activates the target. The concepts of ultrasensitivity, allostery and cooperativity are important not only in understanding the logic of cellular regulation, but also with regard to disease pathology and drug discovery [15].
Classical mathematical models of allostery and cooperative ligand binding, such as the MWC model, were based on observations of cooperativity between symmetric subunits of oligomeric proteins, such as hemoglobin (a tetramer), threonine deaminase (also a tetramer) and aspartate transcarbamylase (a hexamer) [16]. Given that the molecules under study consisted of multiple identical or very similar subunits, it made sense to treat all binding sites as identical. More recently, however, the concept of allostery has been expanded to include monomeric proteins, where binding of a ligand at one site can result in modulation of function or binding at a (perhaps) distant site in the same polypeptide chain [17]. For instance, binding or modification events occurring in an intrinsically-disordered segment of a protein can promote its folding, and this can be communicated to an adjacent segment, with the net effect that a coupled folding-and-binding event or PTM in one region of the protein influences subsequent interactions or modifications at a distant site(s) within the same monomer [18]. Yet another example is hetero-oligomers that display cooperativity such as the ATPase rings in the proteasome and CCT chaperonin complex [19][20][21][22]. In such cases, there is no reason to expect that binding/modification sites will be identical, or that they will make identical contributions to the underlying conformational change once bound/modified.
In the current paper, we set out to explore multisite systems in which the modification of some sites may have a stronger effect on the induced conformational change than the modification of other sites. To do this, we generalize the classical MWC system and assign different parameters to different sites. We aim to determine what combinations of parameters lead to a high level of ultrasensitivity. Each site i is assigned an activation parameter c i , generalizing the parameter c in the original formulation of the MWC model. Small values of the parameter c i correspond to a strong ability for the i-th site to activate the protein. One can also associate to each site a corresponding conformational free energy contribution, that is, the difference in the Gibbs energy function associated to the site i, ΔG i = rt ln(c i ). In other words, rt ln(c i ) is the site-specific free energy contribution to tense-to-relaxed flipping from ligand binding at site i. Notice that the conformational free energy contribution ΔG i is negative number when the activation parameter c i is less than one, and becomes more negative as c i approaches 0. Also note that a large negative ΔG i (and hence a small c i ) corresponds to a strong conformational free energy contribution, which will promote flipping to the active state. In contrast, if c i > 1, then ΔG i will be positive, meaning that the modification does not promote flipping to the relaxed/ active state but instead makes it more likely that the target will stay in the tense state.
Our main results can be summarized as follows. First, making the conformational free energy contribution associated with ligand binding to a single site i more favorable (that is, making this free energy change more negative, which is equivalent to making the activation parameter c i smaller) has a strong tendency to increase the ultrasensitivity of the system, as measured by its associated Hill coefficient H. This effect is not guaranteed as there are some exceptions, especially for low values of the number of sites n, but it holds in most circumstances and under several orders of magnitude for the parameters in the system.
Second, for a fixed number of sites n, each of which is at least moderately active, one can calculate the average of the activation parameters and get a good approximation of the Hill coefficient of the system by assuming that all sites have this average activation parameter. That is, the Hill coefficient is approximately independent on the variability of the parameters c i , only on their mean value.
Third, we find that when the cost of site maintenance is taken into account, one can obtain an optimal ultrasensitive behavior by focusing on a subset of the sites. The strategy is to have a subset of the sites be equally active, and all other sites have a low or negligible conformational free energy contribution. This prediction has been indeed observed in a number of experimental systems, where only a subset of the sites have the ability to activate the protein. In addition, we demonstrate that there are diminishing marginal ultrasensitivity increases in response to conformational free energy contribution improvements, which allows us to predict a maximal effective conformational free energy contribution per site, on the order of -2 to -4 kcal/mol. This prediction follows from first principles of the mathematical model, and it is surprisingly consistent with experimental data for a typical protein phosphorylation site [23][24][25][26].
For completeness, the last sections contain a study of ultrasensitive behavior in a non-allosteric multisite model where all sites are independent from each other, applicable in some cases where the MWC allosteric assumptions are not satisfied. It was found that this system has a more complex relation between activation parameters and ultrasensitivity, which was explored both through computations and mathematical analysis.

Generalized MWC dose response
In this section, we carry out a generalization of the MWC model to account for different activation parameters at distinct sites. See Enciso and Ryerson [8], where a similar generalization was carried out for protein modification efficiencies. Consider a target molecule with n sites in modified-form I, in one of two states, relaxed (R I ) or tense (T I ), where I 2 {0, 1} n is a binary vector representing the modified-form of the target. In the case of protein activation models, relaxed and tense states correspond to different levels of activity. We assume that all modifiedforms in the relaxed (R) conformation are active and have the same activity, whereas all modified-forms in the tense (T) conformation are inactive, and have the same (low) activity. Under MWC assumptions, the relaxed state has a higher affinity to the ligand than the tense state, this is assumed here for the states R I and T I . The unmodified state, I ¼0 ¼ ð0; 0; � � � ; 0Þ, and the fully-modified state, I ¼1 ¼ ð1; 1; � � � ; 1Þ, are the two extreme modified-forms, and a total of 2 n modified-forms are possible. The modification of site i on the target will result in modifiedform J, where J = I [ {i}. In other words J is the modified-form consisting of adding one more modification at site i to the modified-form I. For example, in Fig 1a, a two-site target can be in the relaxed state with no modifications R (0,0) and be reversibly modified to R (1,0) or R (0,1) and subsequently to R (1,1) . A target in a relaxed state can also flip to the tense state in that form. For instance, R (1,1) can flip to the T (1,1) state. Similarly, the tense target in modified-form T (1,0) can be reversibly modified to T (1,1) . We call u the kinase concentration in the case of multisite phosphorylation and u denotes the ligand concentration in the case of ligand binding.
The general system can be described by the chemical reaction network in Fig 1b. The parameter α i is the microscopic association constant for ligand binding at site i. Note that α i is an affinity, i.e. not a dissociation constant but the inverse thereof. If instead of ligand binding the protein is modified by post-translational modification (phosphorylation, acetylation, etc.), α i represents the modification efficiency of site i, which, for example, will be determined by the relative suitabilty of the site to be phosphorylated by a given kinase and dephosphorylated by a given phosphatase. The parameter L is the equilibrium constant between R0 and T0. L is typically assumed to be greater than 1, so that, in the absence of modification, the tense/inactive state is favored over the relaxed/active state. Indeed, L > 9 is required for the target to be less than 10% active in the absence of modification. This network has the property of detailed balance, as the product of the equilibrium constants around any closed cycle of states is 1; this is the same as saying that the net free energy change around any closed cycle of states is 0. For this reason each forward and reverse reaction pair is in equilibrium. https://doi.org/10.1371/journal.pcbi.1007966.g001 We use the notation c I ¼ c i . Notice that Lc I is the equilibrium constant between the relaxed (active) state R I and the tense (inactive) state T I . In this sense, one can think of c i as the contribution to this equilibrium constant made by each individual site i, and ΔG i as the free energy differential between active and inactive protein contributed by modification at site i.
The statistical weights for each state possible when n = 2 is listed in Fig 1c. The probability of a state, say R (1,0) is defined as the ratio between the statistical weight and the partition function Z. In the n = 2 case, Z = 1 + α 1 + α 2 + α 1 α 2 + L + Lα 1 c 1 + Lα 2 c 2 + Lα 1 α 2 c 1 c 2 . S1 Fig  shows a table of the statistical weights of each modification state when n = 3 and the associated We will use u to represent the concentration of ligand or the concentration of modifying enzyme. If the MWC target molecule is a receptor that is regulated by ligand binding, then u is the concentration of ligand. If the MWC target molecule is instead regulated by post-translational modification, then u is the concentration/activity of the modifying enzyme. In the latter case, we assume that the modifying enzyme is in steady state with a corresponding demodifying enzyme (e.g. a kinase-phosphatase system), and that both enzymes are far from saturation. Under these assumptions, the dose response relating kinase concentration u to the fraction of the target sites that are modified is the same as the dose response relating ligand concentration to the fraction of the target sites that are bound. See [6,28] for further details. Following a similar analysis to that of Enciso & Ryerson [8], since the system is in detailed balance, for every index I, Solving for R I and T I , we can relate R I to R0 (relaxed protein with no modifications) by induction as: ðua i þ 1Þ and similarly; Here, r i ðcÞ ¼ X jIj¼i c I is the symmetric polynomial in i with entries c [8]. For example, consider n = 2 with c = (c 1 , c 2 ). Here, ρ 2 (c) = c 1 c 2 , ρ 1 (c) = c 1 + c 2 , and ρ 0 (c) = 1. At various points we are able to rewrite a sum into a product, using the principle that if x i is a constant for i = 1, 2, � � �, n, then ðx j þ 1Þ. For general n, the above allows us to write ðuc j a j þ 1Þ; The response of this system is given by the total concentration of relaxed protein, regardless of its level of modifications. That is, : . Notice from this functional form when any c i is equal to 1, it simply multiplies the dose response by one and becomes the same as a system with n − 1 sites. This is a nontrivial comment which is not obvious from the system otherwise, but it is biologically intuitive. If a target molecule has weak sites, they only contribute weakly or not at all to increase the Hill coefficient. For fixed parameter values c and α, we define the maximal response f 1 ðcÞ ¼ lim u!1 f ðu; c; aÞ. A simple calculation shows that f 1 ðcÞ ¼ 1 1þLc 1 c 2 ���c n and depends only on c 1 , c 2 , � � �, c n and is independent of α. This maximal output value will allow us to normalize response curves across different parameter values in the sections below.
Since the effect of varying the modification parameters α i was extensively described in Enciso and Ryerson [8], in the majority of the discussion here, we will assume that the α i are equal to each other, and in fact we can set α i = 1. To see this, one can re-scale u by defining � u ¼ u� a. The new dose response has the same Hill coefficient as the old system, however the new system satisfies α i = 1 for all i. This helps to better understand the effect of individual activation parameters.

Computational results on MWC ultrasensitivity
Recall from the previous section that f(u, c, α) represents the dose response for the generalized MWC system with parameters c i and α i , for i = 1, 2, � � �, n, where n is the number of sites, and u is the ligand/enzyme concentration.
In this section we carry out a computational analysis of the dose response and its Hill coefficient, where the parameters c i are sampled logarithmically. More specifically, log(c i ) is chosen with uniform distribution between [10 −4 , 0.9]. The parameter L � 1 was fixed, and the parameters α i were chosen to be identical to each other, a i ¼ � a, here the value of � a does not affect the Hill coefficient.
We calculated the Hill coefficient H by solving for EC 90 and EC 10 with a standard numerical solver. Here, we solved for u such that f(u, c, α) − βf 1 (c) = 0 for both β = 10% and 90%. With both EC 10 and EC 90 , we can calculate H as derived in [27]. H > 1 implies the dose response curve is ultrasensitive, while H = 1 implies there is no ultrasensitivity, and H < 1 shows negative ultrasensitivity. One can also think of H > 1 showing that the dose response has a good switch [28]; the larger the value of H the more ultrasensitive the dose response curve.

PLOS COMPUTATIONAL BIOLOGY
Effect of magnitude and variability of energy of activation coefficient tends to increase with the number of sites. In Fig 2b-2d, c 1 and c 2 were increased from 10 −4 to 0.9 and each c i = 0.01 for i � 3, α i = 1, and L = 1000. In these figures, H decreases for increasing c 1 only, suggesting that H increases with increasing conformational free energy contribution (recall that larger values of c correspond to lower activation contributions). Note that for the n = 2 case, there are cases for large values of c 1 where the Hill number is undefined. S2 Fig shows similar data on a linear scale.
In S2 Fig panel e, we show a Monte Carlo approach to study whether H is always a decreasing function of c i . By symmetry, we take any individual c i parameter to be c 1 , without loss of generality. To find the proportion of cases where H decreases on c 1 , we find a numerical approximation to H c 1 ðcÞ as follows: for small Δx to determine if @HðcÞ @c 1  In Fig 3b, for the same parameter values c i and α i = 1 we plot H against the total conformational free energy contribution ΔG tot defined as where r is the gas constant and t is the temperature (traditionally R and T but labeled as r and t, respectively, to maintain consistent notation and not be mistaken for the MWC tense and relaxed states). For fixed n, as the total conformational free energy contribution increases, ultrasensitivity generally tends to increase, and after some threshold, it tends to level off. To increase ultrasensitivity at that point, a target/receptor cannot profitably utilize more conformational free energy contribution, but instead must evolve more sites. As a concrete example, let us consider an MWC molecule with two sites under selective pressure to increase its ultrasensitivity. Changes to the microscopic modification affinities/efficiencies (i.e., the ΔG i 's) will either decrease ultrasensitivity (if the changes are unbalanced), or at best leave ultrasensitivity unaltered (if the changes are balanced) [8]. Thus, the only viable options to increase ultrasensitivity are to (a) evolve another site, or (b) strengthen the conformational free energy contribution of the existing sites. At first, significant increases to ultrasensitivity can result from the second option. A mutation that strengthens the conformational free energy contribution of one of the sites will move the molecule up and to the right in the cloud of points for n = 2 shown in Fig 3b, with the largest jump coming from strengthening the weakest site. As additional mutations of this type arise and become fixed by natural selection, the molecule will move to the top right of the cloud; here the conformational free energies will be approximately balanced and have magnitudes of approximately -2 to -4 kcal/mol. At this point, substantial improvement to ultrasensitivity (i.e., an increase of the Hill number by greater than 0.5 units) can only arise if the molecule evolves an additional site.
To view this more clearly, consider Fig 3c, which shows ultrasensitivity for increasing values of total conformational free energy contribution where parameter values have been set to c i ¼ � c and a i ¼ � a ¼ 1 and fixed n and L. In other words, when each c i is the same, meaning the conformational free energy contribution is the same in all sites, we can see that ultrasensitivity generally increases and eventually levels off as the conformational free energy contribution increases. Fig 3c also helps to make a prediction of the energy that each site optimally contributes, given a total conformational free energy contribution.
Calculating the "knee" of the curve provides a rough estimate of the total conformational free energy contribution required before the ultrasensitivity begins to level off. The knee of a saturating curve is a mathematical definition that captures the point at which the curve is reaching saturation. It is defined in our context as max a � x � b|Y n (c) − ℓ(x)|, where Y n (x) is the Hill coefficient curve in Fig 3c for n sites, a, b are the lowest and highest total energy values among the data points for n sites respectively, and ℓ(x) is a secant line joining (a, Y n (a)) and (b, Y n (b). For a more detailed explanation, refer to Figure 2b in Ref [29].
The approximated knee of each curve in Fig 3c was found and is depicted with a diamond and listed in Table 1. Notice that the energy for saturation increases roughly linearly with the number of sites. In each case, the amount of energy per site is approximately -2 to -4 kcal/mol. This analysis is consistent with some previous experimental findings [23,24]. The approximated H for Ste5 from [24] with n = 8 phosphorylation sites is indicated with an asterisk in Fig  3b and 3c. We derive the value -1.6 kcal/mol per site in this system, which is equivalent to a 10-fold affinity increase per site as approximated in [24]. Not only does this data point approximately lie close to the curve for n = 8, but in fact it lies close to the knee of the curve when all sites contribute an equal amount, as predicted in the above analysis. The marginal effect of an additional kcal of free energy is dependent on only one other parameter, namely L (assuming the α i are roughly equal to each other). If L ranges from 30 to 10,000, the analysis is roughly similar (see S3 Fig panels c-d), and it leads to an energy range of around -2 to -4 kcal/mol per site (see also S1 Table).
There are ways of evaluating ultrasensitivity other than the Goldbeter-Koshland method [27]. In S4 Fig, we measure ultrasensitivity in two additional ways: (1) fitting the dose response curve to the Hill function [9] and (2) Levitzki's n 50 [30]. We see similar results, thus the qualitative results here do not depend on how ultrasensitivity is measured.
In Fig 3d, similar to Fig 3b, we plot H for increasing values of total conformational free energy contribution where now we take into account a maintenance cost for each site, denoted by M c . Such a maintenance cost may arise, for example, if there is rapid turnover of a posttranslational modification, as has been observed for phosphorylation-dephosphorylation of some substrates [31,32]. This type of rapid dynamics in modification-demodification cycles could constitute a non-negligible expenditure of energy for the cell.
The total activation energy including maintenance can be calculated as ΔG tot + M c �n, where ΔG tot is given by (Eq 3) and M c = 4 kcal/mol, which was arbitrarily chosen. In this figure, we assume for simplicity that energy is equally distributed among all sites. Once the cost of maintenance is taken into account, one can see more clearly that for each level of total conformational free energy contribution there is an optimal value of n. For instance, if the total energy is -20 kcal/mol, then the optimal number of sites is n = 3; any fewer sites will not have as high ultrasensitivity, while any larger number of sites requires an excessive amount of maintenance. If there are more than four sites in this system, it is beneficial to eliminate or silence the remaining sites. Similar qualitative results can be seen in S5 Fig with different maintenance cost values.
To summarize, in this section we have shown that (1) increasing the conformational free energy contribution at a single site has a strong tendency to increase the ultrasensitivity of the response, with some exceptions, (2) for fixed n, the ultrasensitivity depends on the mean of the free energies of activation and very little on their variance, and (3) we estimate from first principles an effective energy range of -2 to -4 kcal/mol per site, which is consistent with experimental data.

Generalized independent dose response
The assumption of cooperativity between sites plays a role in the ultrasensitive behavior of the dose response curves. However, if we do not assume cooperativity between sites, will we observe the same effect in the previous section on H? In this section, we use a non-allosteric model and carry out a similar study as for the generalized MWC model. The proposed model has been used elsewhere [8,33] but here it is generalized for the first time to have different activation parameters at different sites.
Consider a target molecule with n modification sites in modified-form I 2 {0, 1} n , where we no longer assume that there is cooperativity between sites. The target can be in one of two states, A I (active) or B I (inactive) and thus gives 2 n possible modified-forms.
The target S in modified-form I can be described by the chemical reaction in Fig 4a, where v i > 1 represents the conformational free energy contribution of the i-th modification site.
Each v i can also be related with the binding energy of the i-th modification site in the MWC model through the formula That is, the larger the value of v i , the larger the free energy. The parameter v i can also be thought of as the inverse of v i ¼ 1 c i of the activation parameter in the MWC model. For notation purposes, v I ¼ To obtain the dose response for total target S T as a function of enzyme concentration u, we use mass action kinetics on the chemical reaction below As in the MWC model, d is the reaction rate constant and is analogous to L. The associated differential equation for the active target is: with conservation of mass equation for the target in modified-form I: We allow this reaction to reach equilibrium by assuming that this activation/deactivation reaction is much faster than protein modification. This is a reasonable assumption in the case of protein phosphorylation. Solving for steady state of (4), In order to calculate the activity level of a target molecule in modified-form I, we defined the function Q I (v), which can be considered to be the fraction of time a protein is active, as To further understand Q I (v), consider the case when n = 2. If both sites are modified, then . In other words, the activity of a protein increases with the amount of modifications. Notice that the activity level will depend on the activation parameters of the specific sites and the overall number of sites modified. We can also determine the concentration of S I , as a function of enzyme concentration u. We can accomplish this by first considering the fraction, p i , that is modified on the i-th site, at steady state. Then, from [8], given k i is the disassociation rate constant of the i-th site, We assume that the modification states of the different sites are independent of each other, an assumption that is in a sense the opposite of cooperativity. In other words, the modification of one site does not influence the modification of another. This allows to calculate the proportion of target in state I as where S T is the total amount of target molecule. The dose response is calculated as follows: This function has a maximal output value f 1 (v), which is found in a similar fashion to the MWC maximal output value by evaluating the limit of f(u, v) as u ! 1. Note that for any I containing a zero (i.e., any modified-form with at least one site un-modified), lim u!1 S I ¼ 0.
Only I ¼1 ¼ ð1; 1; � � � ; 1Þ will contain a non-zero limit for S I , giving f 1 ðvÞ ¼ Q1 . This maximal output value can be used to normalize the dose response curves across different parameter values, similar to the MWC system.

Computational results on independent system ultrasensitivity
For multisite proteins, we can determine the proportion of active target by calculating f(u, v) from the independent system above. In Fig 4b, we plot dose response functions for n = 2, 4, and 8 with v i = 100, k i = 1 and d = 1000.
Similar to the MWC section above, we show how H is affected by the activation parameters of individual sites. In Fig 4c- To determine the effect the variability between parameters v i has on H, we varied parameters v i , measured H and compared to when all parameters v i are equal, similar to the MWC system. In Fig 5a, we sampled a vector with entries from [1, 10 4  The computational and analytical results described in the section below titled "Independent System Mathematical Analysis" suggest that that d > ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p is a biologically reasonable assumption that will give dose response functions where the effect of two modifications is significantly different than the effect of a single modification. Similarly, d < ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p gives dose response functions where the effect of a single modification has a similar effect as two modifications, termed "1+" regime. In this "1+" regime we see H increasing on v. When d ¼ ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p , we have a dose response function where the effect of one modification is approximately 50% of the effect of two modifications, with no ultrasensitivity (H � 1). We can also see that if d is slightly past the 50% of max activation, H can be maximized by increasing the free activation of energy v.
To summarize, in this section we show that (1) ultrasensitivity increases under specific parameter regimes and (2) may depend on the variability between the activation parameters.

MWC system mathematical analysis
In this section, we provide a mathematical analysis of the generalized MWC system showing that H(c, α) is roughly independent of the variation of c. We will show that H is essentially a function of � c and � a. That is, the variability of activation parameters only affects H to the extent that it changes the mean, � c . Consider f(u, c, α) from Eq (1) and define � a ¼ a 1 þ a 2 þ � � � þ a n n ; Da ¼ ða 1 À � a; a 2 À � a; � � � ; a n À � aÞ ; For notation purposes, letĉ ¼ ð� c; � c; � � � � cÞ 2 R n andâ ¼ � a; � a; � � � ; � aÞ 2 R n . Recall that given a C 2 function f such that @f @x ða; bÞ ¼ @f @y ða; bÞ ¼ 0, it holds that f ðx; yÞ ¼ f ða; bÞ þ oðx À a; y À bÞ: We will use this to show that Hðc; aÞ ¼ Hðĉ;âÞ þ oðDc; DaÞ: This formula demonstrates in particular that H essentially does not vary if the mean of c is preserved, as illustrated in Fig 3a. Proposition 1 Hðc;âÞ ¼ Hðĉ;âÞ þ oðDcÞ Proof: For simplicity, assume S T = 1 and assume for now that u andâ are fixed. By the approximation of the geometric mean using the arithmetic mean, we have Taking the n-th power,  One can assume thatâ doesn't change since H is unaffected by increasing or decreasing its value, as explained in the sections above. Also, the above analysis is carried out for u in a neighborhood of EC 10 ðĉ;âÞ and EC 90 ðĉ;âÞ, hence one can assume that u does not vary significantly. Proposition 2 Hðĉ; aÞ ¼ Hðĉ;âÞ þ oðDaÞ. Proof: Similar to Proposition 1, Taking the n-th power,

Independent system mathematical analysis
In this section, the following theorem and proposition provide mathematical analysis for the independent system showing that for n = 2, H increases when d < ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p . Theorem 1 Suppose that f(u, z) > 0 is a saturating C 2 function defined for all u, z > 0, such that f u (u, z) > 0. If the function is strictly increasing on u, then H(z) is increasing for every parameter z > 0. Proof: Let p(z) and q(z) represent the EC 10 kðkÀ uÞ kðkÀ uÞ 10;v 1 þ Q 01;v 1 À �Q 10 À �Q 01 Þ À �kQ 00 kðQ 10 þ Q 01 À 2Q 00 Þ þ uð2Q 11 À Q 10 À Q 01 Þ , so that σ = τ 1 τ 2 . Note that τ 2 is a strictly increasing function on u since t 2 ¼ À C 1 À C 2 u .
In the following text, we show that τ 1 is also strictly increasing on u if and only if d < ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p . To see this, notice that τ 1 is strictly increasing if and only if k C 4 < 1 C 3 , which is equivalent to C 3 k < C 4 . This is equivalent to kð2Q 11 À Q 01 À Q 10 Þ < kðQ 01 þ Q 10 À 2Q 00 Þ 2Q 11 À Q 01 À Q 10 < Q 01 þ Q 10 À 2Q 00 0 < Q 01 þ Q 10 À Q 00 À Q 11 ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p : As long as d < ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p , it follows that τ 1 is an increasing function on u making σ the product of two increasing functions and thus, σ is increasing on u.
Thus, d 2 < v 1 v 2 and so d < ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p . Thus, by Proposition (4), H is increasing as a function of v 1 and v 2 .
The point where d ¼ ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p (the dotted line in Fig 4c) corresponds to the case where the conformational free energy contribution contributed by the singly modified forms is equal to exactly half of that contributed by the doubly modified forms. It can be shown that this situation (which we shall call the 'linear regime') results in a dose response curve with a Hill number of 1. If v 1 and/or v 2 are then increased so that ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p (the region to the right of the dotted line in Fig 4c) becomes greater than d, the singly modified forms now have more than half the conformational free energy contribution of the doubly modified forms, and the Hill number increases. The Hill number will continue to increase until the two singly modified forms, collectively, contribute exactly the same conformational free energy contribution as the doubly modified form. At this point, the system is in the '+ regime', where modification of one, the other, or both sites lead to the same level of activation.
On the other hand, if d > ffi ffi ffi ffi ffi ffi ffi ffi v 1 v 2 p , then the system is closer to the 'both or none regime', where efficient activation only occurs when both sites are modified. Here, increasing v 1 or v 2 reduces the Hill number by pushing the system away from 'both or none' and closer to 'linear'.

Discussion
In a protein with multiple ligand binding sites, the individual sites can differ from each other in two ways: in their microscopic ligand binding affinity, and in the energetic contribution they make, once bound or modified, to functional outcomes such as a ligand-induced conformational change in the bound protein. Likewise, for a protein that is post-translationally modified on multiple sites, the individual sites may have different modification efficiencies, and may also, independently, make differential contributions to downstream functional consequences once modified. For example, in the case of phosphorylation, the amino acid sequence around the target phosphoacceptor residue can substantially influence the efficiency of phosphorylation by the relevant kinase, as well as the efficiency of dephosphorylation by cellular phosphatases [34]. Such tuning of the steady-state level of site modification is biochemically distinct and clearly separable from the effects that the phosphorylation of that site will have on the conformation of the target molecule, its ability to bind other macromolecules, etc. [35][36][37].
Previously, Enciso and Ryerson [8] asked the question "how can the microscopic ligand binding affinities (a.k.a. modification efficiencies) be tuned if the goal is to maximize ultrasensitivity?" Interestingly, they found that ultrasensitivity was maximal when the microscopic affinities were balanced. For instance, for a protein with 4 ligand binding sites, ultrasensitivity was maximized when all 4 sites had exactly the same ligand binding affinity. For a protein with 4 phosphorylation sites, ultrasensitivity was maximized when all 4 sites had the same phosphorylation/dephosphorylation efficiency.
Here we examined how differential energetic contributions of the sites might affect the performance objective of ultrasensitivity. We considered a simple model in which binding/modification promotes a conformational change that flips the modified molecule from an inactive to an active state; this example is readily extended to other known consequences of ligand binding or post-translational modification. We generalized the classic allosteric MWC model to allow for differences in the energetic contributions for any number of distinct sites. We also considered an independent modification model that does not assume allostery or cooperativity among sites. For the generalized MWC system, we found that ultrasensitivity generally increased when the energetic contribution (i.e., the conformational free energy contribution) of one or more of the sites was strengthened. Here, 'strengthened' means that the conformational free energy contribution became more negative; this results in the corresponding activation parameter c becoming smaller. Furthermore, we found that there was no benefit derived from balancing the conformational free energies, nor any penalty for unbalancing them. Our results have implications for understanding the potential trajectories that can be pursued by a protein under selective pressure to increase the ultrasensitivity of its response to modification.
Regarding our finding that decreasing the activation parameter c i of individual sites has a strong tendency to increase the Hill coefficient, this result is analogous to work by Rubin and Changeux [38]. In Fig 2 of that work the authors illustrate computationally that for fixed parameter values of the MWC model, decreasing c leads to an increase in a different version of the Hill coefficient. In our paper we are able to consider individual sites, rather than all sites together, so our result is in a sense a generalization of that shown in [38].
Despite the fact that there is no penalty associated with the conformational free energies being unbalanced, our model nevertheless suggests a factor that may tend to lead to roughly balanced conformational free energies: diminishing returns. Successive, equal-valued improvements of conformational free energy contribution are diminishing with respect to their effect on ultrasensitivity. That is, changes that are of equal magnitude to previous changes increase ultrasensitivity by a smaller amount than the previous changes did. Furthermore, changes to weaker sites increase ultrasensitivity more dramatically than equivalent changes to stronger sites. Eventually, the marginal increase in ultrasensitivity caused by additional improvements to conformational free energy contribution becomes negligible. At this point it can be argued that a zone of effective neutrality has been reached, where the probability of fixation of a new mutation that incrementally improves ultrasensitivity will be essentially indistinguishable from the probability of fixation of a neutral mutation [39]. At this point, substantial improvement to ultrasensitivity can only arise if the molecule evolves an additional site.
There is an additional factor that may further promote the balancing of conformational free energies. Using both computational and mathematical analysis, we showed that when the sites are at least moderately active, the ultrasensitivity is mostly dependent on the mean of the activation parameters and is largely independent from their variance (Fig 3a). Since increasing the variance of the activation parameters tends to be associated with weaker (less negative) total conformational free energy contribution, a prediction is that the sites tend to have roughly equal activation parameters. This can be implemented e.g. by bulk electrostatic mechanisms, which are commonly found experimentally [23,24,40].
These considerations lead to a prediction that conformational free energies will be roughly balanced, with a kcal/mol value roughly equal to the point where ultrasensitivity starts to level out substantially. As shown in Fig 3b and 3c and Table 1 and S1 Table, this "leveling out point" is roughly between -2 to -4 kcal/mol per site, depending on the number of sites n and the level of basal activation (which is determined by the parameter L). This range of -2 to -4 kcal/mol does assume that the efficiency of modification is roughly constant across all sites, but is otherwise surprisingly independent of other parameters. For example, the range found changed very little upon variation of the value of L from 30 to 10,000, therefore covering most biochemically realistic values for this constant. The range of approximately -2 to -4 kcal/mol corresponds to activation coefficient (c) values between approximately 0.05 and 0.001. Such c values are all within the range reported for classic "MWC enzymes" such as threonine deaminise, glucose-6-phosphate deaminase, aspartate transcarbamoylase and glyceraldehyde-3-phosphate dehydrogenase [41][42][43][44]. Moreover, with regard to phosphorylation, the effect of a single phosphate on conformation [25], protein-protein binding [26] or protein-membrane binding [23,24] has been estimated to be about 2 kcal/mol.
We also showed (Fig 3d) that when the number of sites is large, and a hypothetical maintenance cost per site is included (such as might arise from rapid phosphorylation-dephosphorylation cycles [31,32]), an optimal strategy to maximize ultrasensitivity can be to focus on a subset of the sites, and essentially keep the other sites silent. In such cases, evolving another site is not a viable strategy to increase ultrasensitivity, and it can be argued that there is an optimal number of functional sites that will maximize benefit (ultrasensitivity) while containing cost.
This analysis applies for other forms of multisite modification other than phosphorylation such as ligand binding, methylation, acetylation, etc. When the multisite target molecule has a symmetric structure (such as hemoglobin which is a tetramer), one can assume that the conformational free energy contribution is similar across all sites. In this sense the current study is most relevant when the target structure is more heterogeneous, such as in the case of phosphorylation. Although phosphorylation consumes energy and is not thermodynamically closed, the MWC model is still a popular model to describe it [8,45,46]. It is also mathematically more amenable than the non-allosteric, independent model that we also included for completeness.
Work by Kafri et al [19] has previously studied a mathematical model of chaperon-containing TCP-1 protein that has several sites with different ATP binding affinities. This system can provide very interesting parallels with our framework. Their mathematical model shows that when a protein has multiple sites with different ligand affinities, the Hill coefficient can be reduced leading to apparent negative cooperativity. We do observe a similar effect (see eg Fig 3), although that model has important differences such as variability in modification affinity rather than conformational free energy contribution.
The analysis in this manuscript is limited to systems in equilibrium, i.e. the long term response to a constant input. For non-equilibrium systems, and in situations where energy is used, recent work by Estrada et al. [47] shows that one can obtain a larger Hill coefficient. The authors use techniques similar to kinetic proofreading, which can give rise to large response differences given small differences in ligand affinity. A full discussion of non-equilibrium dynamics is however outside of the scope of our work.
Many dose response curves for allosterically-regulated proteins can be well-modeled by the standard MWC model. Our goal in generalizing the MWC model was to explore the qualitative theoretical consequences of allowing the conformational free energies of different sites to vary, and not to make a tool for empirical fitting to data. On this point, however, it should be noted that Stefan et al. [48] have shown how an extended MWC model such as the one developed here can be used in parameter estimation, and experimental methods to measure MWC parameters are constantly improving [49][50][51]. Another useful recent tool is a method developed by Gruber et al. [52] which facilitates the determination of parameters in the MWC model, by finding a theoretical relationship between the Hill coefficient and model parameters.
Supporting information S1 Fig. Statistical weights of MWC modification states for n = 3. Table of each possible modification state in the generalized MWC system above when n = 3 and the corresponding statistical weight of that state. When n = 3, the associated partition function Z = 1 + α 1 + α 2 + α 3 + α 1 α 2 + α 1 α 3 + α 3 α 3 + α 1 α 2 α 3 + L + α 1 c 1 L + α 2 c 2 L + α 3 c 3 L + α 1 α 2 c 1 c 2 L + α 1 α 3 c 1 c 3 L + α 2 α 3 c 2 c 3 L + α 1 α 2 α 3 c 1 c 2 c 3 L.  [30] as H Lev = 4 � EC 50 � f 0 (EC 50 , α, c) where f 0 (EC 50 , c, α) is the derivative of the dose response function evaluated at the EC 50 , the effective enzyme/ligand concentration at which there is a 50% maximal protein response, labeld H Lev . We can consider H Lev as the sensitivity at 50% maximal response. EC 50 was found with the standard MatLab fzero solver and the derivative with diff after normalizing to the f 1 (c).  Table. Ultrasensitivity at knee. Ultrasensitivity as measured by the Goldbeter-Koshland formula described in Eq (2) along with the approximated knee of curves similar to those in Fig  3c for fixed values of L and n. Parameters a i ¼ � a ¼ 1.