Compact Modeling of Allosteric Multisite Proteins: Application to a Cell Size Checkpoint

We explore a framework to model the dose response of allosteric multisite phosphorylation proteins using a single auxiliary variable. This reduction can closely replicate the steady state behavior of detailed multisite systems such as the Monod-Wyman-Changeux allosteric model or rule-based models. Optimal ultrasensitivity is obtained when the activation of an allosteric protein by its individual sites is concerted and redundant. The reduction makes this framework useful for modeling and analyzing biochemical systems in practical applications, where several multisite proteins may interact simultaneously. As an application we analyze a newly discovered checkpoint signaling pathway in budding yeast, which has been proposed to measure cell growth by monitoring signals generated at sites of plasma membrane growth. We show that the known components of this pathway can form a robust hysteretic switch. In particular, this system incorporates a signal proportional to bud growth or size, a mechanism to read the signal, and an all-or-none response triggered only when the signal reaches a threshold indicating that sufficient growth has occurred.


S1.1: Analysis of the Zds1 -PP2A interaction subsystem
Recall that Z, P are the active monomer concentrations of Zds1, PP2A respectively, Z , P the total monomer concentrations, and p , z the fraction of modified sites. Assume the From the MF framework, estimate Z = Zh 1 (z), P = Ph 2 ( p), where h i (x) = e γ i x / (δ i + e γ i x ) . The Zds1/PP2A dimer is assumed to be active exactly when both monomers are active, so if D is the active form concentration D D = Z Z P P = h 1 (z)h 2 ( p), that is, D = Dh 1 (z)h 2 ( p).
The differential equations for this system are: Notice that the dynamics of z depends on D , the active dimer. The dynamics of p depends on the concentration of the protein C , the active Rho1/Pkc1 dimer, which is treated for now as a given constant. One can eliminate the variables z in , p in and replace them with z, p using the conservation equations z in + z = 1, p in + p = 1. Similarly, one can use the mass conservation of Zds1 and PP2A, Z + D = S Z , P + D = S P to eliminate Z and P from the system. Here S Z ,S P are the total amounts of Zds1 and PP2A including monomer and dimer forms. Finally, the active dimer D can be written in terms of the other variables as described above. The new equations are Now, over time the first and third equation of this system converge towards a constant independent of the other variables. Solving for p at steady state leads to . The quadratic equation for D at steady state has two solutions, however one of them is larger than S Z , S P and can be discarded. Thus over time D converges towards a steady state value D 0 independent of C .
Putting all together, the system is reduced to a single differential equation for z . As stated in the main text, the solutions of the system converge towards the steady state whenever z does, and one can obtain bistability depending on parameter values. The differential equation for z can be rewritten as Both terms in the right hand side are described in Figure 3B as separate functions, and the steady states correspond to the intersection of the two graphs. By comparing both functions on the left and right hand sides of each steady state, one can determine that the first and third steady states are stable, and that the second steady state is unstable. This is a general property of the system whenever there are three steady states, since the sigmoidal shape of h 1 (z) is preserved under other parameter values, and therefore when multiplied by the constant h 2 ( p) and the function (1− z) , the qualitative shape of the graph is preserved. Therefore one can conclude that the Zds1/PP2A interaction system is bistable whenever there are three steady states.
The canonical output of the core model is the active PP2A/Zds1 dimer concentration D . This value can be calculated from the equation above, but more easily from the steady state values of z using the equation . This is the equation used to plot the graphs of C versus D in Figure 3C.

S1.2: Analysis of the Rho1 -Pkc1 interaction subsystem
This system involves the proteins Rho1, Pkc1, denoted R , K respectively in their active forms. Total Rho1 and Pkc1 monomer concentration, regardless of whether active or inactive, are R, K , while k denotes the fraction of modified Pkc1 sites. The given chemical reactions are The MF approximation is used to calculate A few conservation of mass equations hold: Using mass action kinetics one can derive a system of differential equations describing the behavior of this system. In the following, the variables k in and K are eliminated and replaced by 1− k and S K − C , respectively: A steady state analysis of this system is now carried out. Even though the system contains multiple variables, it is still possible to reduce all equations to a single 1D equation at steady state. After setting all left hand sides to zero, notice that the equation for k only depends on C (and the fixed input parameter D ). One can solve it to find where we introduce the notation Q i = q − i / q i for every i .
Next, the steady state equations for R in and R depend linearly on the variables R in , R .
The determinant of the linear system is the denominator of the equation, and it is larger than zero (S K − C = K ≥ 0 ) so that the solution is always unique. For convenience, call Y the denominator of this expression. Replacing into the equation for C , we have Y dC dt = q 4 (S K − C)[q 6 q −4 C + q −9 q −4 C + q 6 q 9 v] − Yq −4 C. By multiplying out on the right hand side, the first two terms of each summand cancel out, so that Y dC dt = q 4 q 6 q 9 (S K − C)v − q −4 q −6 q −9 C + q −4 q −9 C(q 6 + q −9 ). Therefore Y q 4 q 6 q 9 By multiplying on both sides by h 3 (k) and using the identity C = Ch 3 (k) , one obtains Finally, replacing k by its expression in terms of C , at steady state In order to determine the stability of the steady states, notice that h 3 (k)Y / q 4 q 6 q 9 > 0 , even though this expression might depend on k and C . By comparing the sign of the two summands in the right hand side of the expression for dC / dt , one can determine when the steady states are stable or unstable. In the example in Figure 3D, once again one can see that out of the three steady states, two are stable and one is unstable. S1.3: Rule--Based Derivation of MWC We begin by introducing in detail the assumptions that lead to the MWC model in the context of multisite phosphorylation. Recall that a substrate P is assumed to be phosphorylated nonsequentially by a kinase E , that the phosphorylation is cooperative, and that the number of sites is large enough that each individual phosphorylation has a small effect on substrate activation. The substrate may be either active or inactive at any given time, due to structural transitions or the binding of the substrate to another molecule. Based on these assumptions, we start with a basic model containing a large number of variables and reduce it to the MWC model. This framework also relates our modeling framework to the approach known as rule--based modeling, where a series of chemical reactions is defined using a streamlined algorithm, and high--powered computing is used to handle the resulting large number of variables.
Suppose that J ∈{0,1} n is any vector of n zeros and ones indicating the phosphorylation state of the substrate, also known as its phosphoform. A J and I J indicate the concentration of the active and inactive protein phosphoforms in state J , respectively. For any J , K ∈{0,1} n , we say that J  K if J i ≤ K i for every i and | K |=| J | +1, that is the K --phosphoform has exactly one more phosphorylation than the J --phosphoform. Define the chemical reactions displayed in Figure S1A, where J is any fixed index and K is such that K  J . That is, every A J can transition into I J and back, and A J can also become phosphorylated at any of its remaining locations at a rate proportional to E . Every additional phosphorylation makes the transition of active to inactive substrate slightly more difficult, since the transition rate is multiplied by ε |J| , where ε < 1 . Since each phosphorylation is assumed to have a small effect, we can assume that ε is close to 1.
Often there are so many reactions in such an approach that it is not possible to carry out a mathematical analysis, but we can do it in this case. Putting together all reactions, we have for every index J : The idea is to cluster together many variables A J by carrying out a change of variables. Define Then The rate equations for A 0 , A n , and I i for i = 0…n can be calculated in the same way. The full system has the following equations: This is exactly the system of equations for A i , I i in the diagram shown in Figure 2B. In this way, MWC is a sequential shorthand model for a nonsequential system with 2 × 2 n variables. The same framework can be used if the sites are modified through any covalent enzymatic reaction such as acetylation or methylation. They could also be non--covalent ligand binding sites, such as a receptor complex, with a background ligand E in large enough amount that its concentration is not affected by its binding to the receptor. Recall that to calculate the activation function h(x) , in the main text we assumed an approximate balance at steady state for the following reaction of activation and inactivation after i phosphorylations: After calculating the fraction of active sites with this number of phosphorylations, this shows that where γ = −nlnε and δ = L 1 / L 2 . The balance between A i and I i is guaranteed for the parameters as defined above since the system MWC satisfies the property of detailed balance [1]. However it can also be approximately satisfied when detailed balance fails, especially when the transition rates L 1 , L 2 are large compared with other parameters. Detailed balance is usually satisfied for systems that conserve energy such as multisite ligand binding. Systems that require energy to function, such as phosphorylation which requires ATP, are known as "dissipative" and don't necessarily satisfy this property, in which case the transition rates may play a more important role when compared with other rates in the system. In the case of the cell cycle checkpoint proteins, the standing assumption is that their activity is determined by internal transitions or by binding to a regulatory protein, which would need to be compared with e.g. the rate of the different enzymatic reactions.