Abstract
Background
Translation efficiency of certain mRNAs can be regulated through a cytoplasmic polyadenylation process at the preinitiation phase. A translational regulator controls the polyadenylation process and this regulation depends on its posttranslational modifications e.g., phosphorylation. The cytoplasmic polyadenylation binding protein (CPEB1) is one such translational regulator, which regulates the translation of some mRNAs by binding to the cytoplasmic polyadenylation element (CPE). The cytoplasmic polyadenylation process can be turned on or off by the phosphorylation or dephosphorylation state of CPEB1. A specific example could be the regulation of Calcium/Calmodulindependent protein kinase II (αCaMKII) translation through the phosphorylation/dephosphorylation cycle of CPEB1.
Result
Here, we show that CPEB1 mediated polyadenylation of αCaMKII mRNA can result in a bistable switching mechanism. The switch for regulating the polyadenylation is based on a two state model of αCaMKII and its interaction with CPEB1. Based on elementary biochemical kinetics a high dimensional system of nonlinear ordinary differential equations can describe the dynamic characteristics of the polyadenylation loop. Here, we simplified this highdimensional system into approximate lower dimension system that can provide the understanding of dynamics and fixed points of original system. These simplified equations can be used to develop analytical bifurcation diagrams without the use of complex numerical tracking algorithm, and can further give us intuition about the parameter dependence of bistability in this system.
Conclusion
This study provides a systematic method to simplify, approximate and analyze a translation/activation based positive feedback loop. This work shows how to extract low dimensional systems that can be used to obtain analytical solutions for the fixed points of the system and to describe the dynamics of the system. The methods used here have general applicability to the formulation and analysis of many molecular networks.
Background
Cellular signaling pathways that can operate in a switch like manner are called bistable systems [1,2]. A bistable system has the ability to switch between two distinct stable steady states and such a system cannot rest in any intermediate state [3,4]. In response, to an external stimulus a bistable system can move from one state to another. If this switching is permanent then such a system is called irreversible otherwise it is a reversible switch [5]. The bistability in a signaling network is typically due to a positive feedback loop or double negative feedback loop [3]. However, the presence of a positive or a double negative feedback loop does not guarantee bistability [4]. In addition, to these feedback loops a biological network must have nonlinear interactions to exhibit a bistable behavior. Previous, experimental work has described several examples of naturally occurring bistable system [618]. Still bistability is not considered to be a unifying theme of cellular signaling networks and more experimental work is needed to establish bistability as one of the general mechanism of cell signaling [4]. Typically, bistable biological systems were described either at the level of gene expression due to the regulation of gene expression by transcription factors, or at the level of posttranslational modifications e.g., activationdeactivation cycle due to phosphorylation. Here, we present a model of bistability that can arise from the control of gene expression at the level of translation of new proteins
Cytoplasmic polyadenylation regulates the translation efficiency of certain mRNAs through modulating the length of 3' poly (A) tail [19,20]. Polyadenylation can be regulated through a short nucleotide sequence known as cytoplasmic polyadenylation element (CPE) located in their 3' UTR. A CPE binding protein (CPEB1), represses the translation through its dual interactions with CPEs and other mRNA binding proteins. Phosphorylation of CPEB1 changes its interactions with these other proteins, promotes polyadenylation and increases the length of the poly (A) tail [19,20]. The mRNAs with a longer poly (A) tail are more likely to be translated compared to mRNAs with a shorter poly (A) tail [21,22]. The exact mechanism through which the longer poly (A) tail enhances the translation efficiency is not clear. However, it is believed that longer poly(A) tail leads to a circular mRNA which enhances the translation efficiency through recycling the translation machinery on the mRNA frame thus increasing the possibility of mRNA translation through initiation[23,24].
The translation of a highly abundant brain protein αCaMKII (2% of total brain protein is αCaMKII) is regulated through activity induced polyadenylation. The αCaMKIImRNA contains two CPE elements in its 3'UTR and its translation can be regulated through phosphorylation of CPEB1 [25,26]. Recent studies have shown that αCaMKII can phosphorylate CPEB1 and therefore, possibly modulate its own translation through a positive feedback loop. Here, we examine the hypothesis that the positive feedback loop between αCaMKII and CPEB1 forms a bistable switch which regulates the translation of αCaMKII. The aim of this paper is to obtain analytical expressions for the bistability of the CPEB1 αCaMKII molecular pair as a generic example for such systems.
In this paper, we analyze the mathematical properties of this molecular loop. The characteristics of this loop are analyzed by evaluating the dynamics and directly locating the fixed points. Using the elementary biochemical kinetics we develop a molecular model of selfsustained polyadenylation based translation loop. This simple molecular model is represented by six differential equations. We simplify this model through introducing an approximation and algebraic manipulations. These, systematic simplifications result in a three dimensional differential equation based model. Further approximations can reduce this to a single dynamical equation. Based on our approximate equation we also developed an approximate analytical method to directly locate the fixed points of this system. We compare these analytical results to the numerical bifurcation diagrams obtained through numerical tracking of the complete system of equations and show their correspondence. Our results demonstrate that such a positive feedback loop which involves the control of translation through polyadenylation can indeed be bistable over a wide range of parameters. This simplified model, though motivated by the αCaMKIICPEB1 loop, could be seen as a generic model for such a positive feedback loops which involves translation, and degradation of proteins. Such feedback loops do not conserve the quantity of these proteins and are therefore qualitatively different than most post translational feedback models [3,4]. Since, we use a simplified model we can obtain approximate analytical results, which provide us with intuition about how such feedback systems operate.
Method
A. Complete Model Equations
The following set of reactions is used to describe the interactions between the αCaMKII and CPEB1 molecule. These biochemical reactions are based on elementary MichalisMenten type kinetics. The dynamical variable X represents αCaMKII and Y represents the CPEB1. The P subscript represents the phosphorylated form and an A as superscript represents the active form.
From above reactions following differential equations can be deduced.
B. Model analysis and reduction
In the following equations, the dynamical variable X represents αCaMKII and Y represents the CPEB1. The P subscript represents the phosphorylated and active form. By using the pseudosteady state assumptions the differential equations representing the complexes [C_{1}C_{3}] can be eliminated:
The unphosphorylated CPEB1 (Y) is related to phosphorylated CPEB1 (Y_{P}).
Where, Y_{T }is the total amount of CPEB1. The αCaMKII molecules are either in free or in bound form, therefore, the total concentration of αCaMKII is given by following equation.
The differential equation representing the phosphorylated CPEB1 (Y_{P}) from R_{1}R_{6}
By requiring a steady state and substituting we get
Also from (B1,B2, B3, B4, B5 and B7) we can obtain the value of X.
Where, N and M are constants defined as
As described in the result section the equation 2 is as follows:
Substituting C3 from (B3) in [2]
Placing Y_{P }from (B7) into (B9)
Where, P_{7 }= P*k_{7 }and P_{4 }= P*k_{4}. Defining, U = (Ca^{+2})_{4}CaM, and c = b*Y_{T }or
We obtain
The equation 1 describes the rate of change of total concentration of αCaMKII (X_{T}) is obtained from the balance between the new synthesis of αCaMKII and its degradation. This equation is further transformed to the approximate solution as shown by equation 4. Similarly, from R_{1}R_{6 }an equation representing the dynamics of phosphrylated αCaMKII can be constructed.
The equation B12 is further simplified in only two variables i.e., X_{T }and X_{P }by placing the and putting the values of X, C_{1}, C_{2}, Y, and Y_{P }from B1 to B8 and further simplifying we get the following expression.
The I_{1},I_{2},I_{3 }and I_{4 }are defined as follows:
The expression in B14 to B17 are defined as k_{7 }= k_{7 }*P, k_{4 }= k_{4 }*P,a_{1 }= (k_{3 }*k_{1 })./(k_{2 }+ k_{3 }) a_{2 }= k_{10 }*U, D_{I }= k_{4 }+ k_{1010 }, d_{1 }= (k_{6 }.*k_{5})./(k_{55 }+ k_{6})
The expression of B13 gives us a function where X_{T }= f(X_{P}), It is almost impossible to invert this function to get X_{P }= f^{1}(X_{T})so we numerically approximated it through function fitting.
Equation B12 with the aid of equation B1B 13 can be transformed to following polynomial in terms of X_{P}.
Where, the coefficients z_{12},z_{13 }and z_{14 }are defined by following expressions (B19B32). The equation (B18) is obtained through further simplifying (B13) such that a fourth order polynomial is obtained in terms of X_{P}. This 4'th order polynomial has an analytical solution because it has no zero order term and therefore has one solution X_{P }= 0, and the other solutions are the solutions of a third order polynomial.
Results
Analysis of the model of αCaMKII synthesis through polyadenylation
Our simplified model of a selfsustained polyadenylation of αCaMKIImRNA (Figure 1) is based on biochemical interactions between a plasticity related kinase αCaMKII and its translational regulator CPEB1 through a positive feedback. This model of polyadenylation based translation of αCaMKIImRNA (Figure 1) is composed of two molecular components which interact through a closed loop. 1) The αCaMKII protein which can be in two states: inactive, and phosphorylated, active. 2) The CPEB1 a translational regulator to regulate the polyadenylation, can be either in the phosphorylated or unphosphorylated state. The phosphorylated CPEB1 promotes the translation at preinitiation phase through polyadenylation. Here, in this simple model we assume that CPEB1 is phosphorylated only by active and phosphorylated αCaMKII. In this molecular scheme the αCaMKII protein is removed at a certain degradation rate [27]. The synthesis of new αCaMKII protein regulated by polyadenylation provides the necessary compensation for the amount removed due to degradation. Our model shows how the concentration of αCaMKII can be maintained at multiple levels despite the synthesis of new molecules and removal due to protein degradation.
Figure 1. The simplified model of CPEB1 mediated polyadenylation of αCaMKII through a selfsustaining αCaMKIICPEB 1 molecular loop. The αCaMKII, molecule can be inactive, or active and phosphorylated state, while CPEB1 can be in unphosphorylated and phosphorylated states. The αCaMKII molecule in both these states degrades with a certain degradation rate. The active and phosphorylated CaMKII phosphorylate the CPEB1 molecule, which through polyadenylation of αCaMKII mRNA generates the new CaMKII molecule. Here, the synthesis of new αCaMKII molecule is represented by a translation step (T) which, describe the multistep translation process through a single biochemical reaction R6 as described in appendix.
This system is described by a set of differential equation (Method A). We use notation in which CaMKII is denoted as X, and phosphorylated αCaMKII as X_{P}.
Similarly CPEB is denoted as Y and phosphorylated CPEB1 as Y_{P}. On purpose we have chosen a simplified version of these components, when compared to our previous work [28] in order to gain intuition into the behavior of the system while keeping the bistability. We can easily show that for appropriate parameters this system is bistable. However, gaining an intuitive understanding of this bistability is difficult even in this simplified system, because the dimensionality is still too large.
The twostate assumption of the αCaMKII is a significant simplification compared to more realistic multistate models and was chosen to enable simpler analysis. The synthesis of new proteins in our proposed molecular model is controlled by the phosphorylated CPEB1 molecules. The translation model is also a simplified representation of a complex molecular system. The aim of this paper is to generate a much reduced dynamical model that can nevertheless capture the qualitative behavior of the system, and provide us with a more intuitive understanding of dynamical behavior of the system.
One level of simplification is to set the derivatives of the three complexes (C_{1}C_{3}) to zero (Method B), to obtain a pseudo steady state approximation, similar to that used in the standard michaelismenten approximation. This reduces our system to 3 dynamical equations, and three functional expressions, as described by equations 1, B6, and B12. This approximation simplifies the dynamics, and can be used to more easily find the steady states of this system numerically. We can further approximate this system by a single differential equation. The key simplifying assumption required here is the following equation:
Where, C_{3}, as described by equation B3 in the method B, is the concentration of the phosphorylated form of CPEB1, bound to the translation machinery, or in other words the concentration of active translation machinery available for producing new αCaMKII. Here k_{9 }is the forward rate of generating new αCaMKII. This equation assumes that αCaMKII degrades at the same rates both in free and bound forms. The X_{T }in above equation (Eq. 1) is the concentration of total αCaMKII. In our detailed system of equations (Method A) we do not assume the degradation of complexes. However, this approximation holds if the relative concentration of complexes at steady state is small.
By using a pseudo steady state approximation on some of the faster dynamical variables, we obtain (Method B) that:
Where the parameters 'c', 'P7' and 'a', are defined in Method B. The parameter c is linearly dependent on Y_{T}, and is inversely proportional to degradation rate λ. It shows that amount of total αCaMKII increases with increase in translation (either due to increase in translation rate or concentration of translation machinery "T") and decreases with increase in degradation rate. The K_{m1/2 }of this process is P_{7}/a and depends on phosphatses and dephosphorylation rate of CPEB1.
We call the first function of RHS of equation 3 the synthesis function and the second one a degradation function. We can rewrite equation 3 as
Where, G is a synthesis function, and F a degradation function, and the dynamics are then understood as a balance between synthesis and degradation. If we can derive a function that relates X_{p }to X_{T}, i.e., X_{P }= g (X_{T}), then equation 4 could be converted to equation of a single dynamical variable:
and we would have reduced the 6 coupled ODE's to a single ODE:
Although we have an exact expression for X_{T }as a function of X_{P }[(Method B equation 13)], it is difficult to analytically invert it to obtain an analytical expression for X_{P }= g ( X_{T}). Instead it is easy to invert this numerically using a higher order polynomial fit.
Using this fit.
The fitting performance of this equation is shown in additional file 1, Figure S_{1}.The solid black line in additional file 1, Figure S_{1 }represents the original function {Method B B13}, while the dotted red line represents the fit from equation 7. Thus we reduce the 6 differential equations to a single approximate ODE of the form described in equation 6, where:
Additional file 1. Figure S1 The fitting performance of equation7 with original function B13 (Method B).
Format: PDF Size: 16KB Download file
This file can be viewed with: Adobe Acrobat Reader
Here, the X_{P }is obtained from fitting function 7 and F (X_{T}) is given by following equation.
In order, to gain some intuitive understanding of this system we graphically plot the two terms in equation 6 (Figure 2a). The function G is shown as dotted red line and F as solid blue line. This graphical representation (Figure 2a) provides a very simple intuitive explanation of the behavior of the system. The function G' can be seen as the synthesis function (a source) and the function F as the degradation function (a sink). The intersections between these two curves are the systems fixed points. For low values of X_{T }there is hardly any synthesis, and degradation dominates, so the system converges to the low fixed point at zero. At higher values of X_{T }there is an abrupt rise in the level of protein synthesis, which quickly saturates. The cross between G and F in the quickly rising portion of G is the unstable fixed point, below it degradation dominates and above it synthesis dominates, and the level of X_{T }increases until it reaches the third intersection, which is the upper stable fixed point. Due to the saturation of G, at higher levels of X_{T}, degradation dominates again, resulting in a convergence from above to the upper stable fixed point. Changing the degradation rate will simply change the slope of the F function, and changing synthesis and activation parameters, will quantitatively change the shape of the G' function.
Figure 2. The approximate analytical solution of αCaMKIICPEB1 molecular loop. It is obtained through graphically locating all the steady state solutions of single differential equation 1. This equation describes the rate of change of total concentration of αCaMKII (X_{T}) and its functional relation between the new synthesis and degradation of CaMKII. The steady state form of equation 1 is transformed into two functions F (X_{T}) and G (X_{T}). The function F (X_{T}) (solid blue curve) basically represents the net αCaMKII degradation, whereas the modified function G'(X_{T}) = G [g (X_{T})] (dotted red curve) represent the net generation of αCaMKII. (a) This analytical solution is developed by setting all parameters as described in table 1. Graphs of both functions intersect at three locations which are characterized as steady state solutions of equation 1 and provide the approximate analytical solution. The upper steady state solution is located at X_{T }= 95, while the unstable solution is located at X_{T }= 9.4, and lower steady state solution at X_{T }= 0.0001. (b) The effect of degradation rate on analytical solution. Four different degradation rates are selected λ = 0.00006 s (solid blue line #1, bistable system), λ, = 0.0001 s^{1 }(solid blue line #2 bistable system), λ = 0.0003 s^{1}(solid blue line #3 bistable system) and λ = 0.0006 s^{1 }(solid blue line #4 monostable system). (c) The effect of CPEB1 activation parameter k_{5 }on analytical solution. Four different values of CPEB1 activation are selected k_{5 }= 0.00001μM^{1}.s^{1 }(dotted red line #1, monostable system), k_{5 }= 0.00006μM^{1}.s^{1 }(dotted red line #2 bistable system), k_{5 }= 0.0001μM^{1}.s^{1 }(dotted red line #3 bistable system) and λ = 0.0072 μM^{1}s^{1 }(dotted red line #4 bistable system).
Next we show the impact of changing some of the system parameters on the steady state solution characteristics of equation 6 (Figure 2b and 2c). Changing the degradation rate only affects the degradation curve (Figure 2bsolid blue line). As the degradation rate is increased from slower (Figure 2bsolid blue line 1) to much faster (Figure 2bsolid blue line 4) the system moves from robustly bistable to a monostable system. Thus, at much faster degradation the amount of CaMKII generated is not enough to compensate for the loss due to large degradation rate therefore, system has only lower stable steady state solution (Figure 2bsolid blue line 4) as the CaMKII degradation rate is decreased a balance between new synthesis and degradation is restored and system becomes a bistable switch (Figure 2bsolid blue lines 24). As the degradation rate is decreased, and the system as shifts from monostable to bistable, the lower fixed point does not change, the value of X_{T }at the unstable fixed point is increased, and the value of X_{T }at the upper fixed point first rapidly increases and then plateaus (Figure 2bsolid blue line 2,3,4).
Other parameters affect only the synthesis curve (G'). One such example is parameter k_{5}, which quantifies the CaMKII mediated activation of CPEB1 (Figure 2c). When this activation parameter is set at low value (Figure 2cdotted red line 4) the system has only lower stable steady state solution. However, as the activation rate of CPEB1 is increased, more active CPEB1 is available, which in turn stimulates new protein synthesis thus shifting the synthesis curve in upward direction (Figure 2cdotted red line 3,2,1) and shifting the system from monostable to robustly bistable state.
If we are only interested in the steady states of this system we can set all the derivatives to zero and obtain the following forth order polynomial to describe the steady states of the system:
Where, the coefficients z_{12}z_{14 }are defined in method B. The derivation leading to this result is explained in method B. This equation is more precise than equation 6 because here we avoided the need to use the polynomial fit for the inversion. However, this equation can only account for the fixed points, and is unable to account for system dynamics. This 4'th order polynomial shows that this system has 4 steady states, one is always at zero because there is no zero order term. The signs of the coefficients determine how many real solutions this polynomial could have [29] and biologically we are only interested in positive real solutions.
Comparison of the different approximations to the full model
Our aim here is to describe the behavior of the synthesisdegradation loop (Figure 1) and the various approximations used to simplify it. We can integrate the dynamics either using the full system of ODE's (Method A, A1A6) or by using the reduced three dimensional systems (Method B, equations 6, B6, B12), or the approximate one dimensional system (equation 6). The fixed points, and their parameter dependence can be found either by numerical bifurcation analysis of the full system (Method A, A1A6), the reduced three dimensional system (6, B6, B12), the approximate 1D system (equation 6) or by the 4'th order polynomial [10]. The first three methods can be used for obtaining the system dynamics and the last method only for the fixed points.
In order, to further compare the three dynamical equation methods we first analyzed the dynamical properties of this system and then the steady state solution characteristics of this loop with respect to certain system parameters. The dynamics of this system (Figure 3) is analyzed with two different set of parameters. The first set of parameters (Condition I) represent a case where the approximation of equation 1 is a good aproximation (Figure 3a), whereas the second set of parameters (Condition II) represents a case where this approximation does not hold (Figure 3b). For condition I the approximation holds since the amount of αCaMKII bound in biochemical complexes is negligible, whereas for condition II a significant portion is trapped in biochemical complexes. We obtain these different conditions by setting the different values of autophosphorylation parameter k_{3}. For condition I we set k_{3 }= 500μM^{1 }s^{1}, whereas for condition II we set this parameter at 0.5 μM^{1 }s^{1}. The values of all the other parameters are identical in both conditions (Table 1). For all these cases a 10 second (Ca^{+2 })_{4}CaM pulse is used to provide the necessary stimulus to drive system from the basal to the upstate.
Table 1. The numerical values of different parameters of bistable molecular loop.
Figure 3. The dynamic characteristics of aCaMKIICPEB 1 molecular loop. These characteristics are developed through three alternative models. First model is based on differential equations A1A6 (Full model, represented by solid red line (up state) and solid blue line (down state)). Second model is based on differential equation 1, B6 and B9 (Reduced model, represented by dotted red line (up state) and dotted blue line (down state)). Third model is based only on a single differential equation (equation 6) (Single differential equation model, represented by dotted black line (up state) and dotted blue line (down state)). System moves from down to upregulated state through a 10 second application of (Ca^{+2 })_{4}CaM pulse. The dynamic simulations are carried out with two different set of parameters. For one set of parameters the dynamic from three alternative models converge to same upper and lower steady states (a). In this case the dynamics converge to same upper steady state at X_{T }= 95 and lower steady state at X_{T }= 0.0 For second selection of parameters the dynamic results from three alternative models does not converge to same upper steady state (b) although the lower steady state remains the same for all three models. Here, the dynamics from full model converge to upper steady state at X_{T }= 360 and lower state at X_{T }= 0.0001 (Figure 3 b, solid red line (up state) and solid blue line (down state)), whereas the dynamics from reduced model (dotted red line (up state) and dotted blue line (down state)) and single differential equation model (dotted black line (up state) and dotted blue line (down state)) converge to upper steady state at X_{T }= 95 and lower state at X_{T }= 0.0001.
For condition I the dynamics from all three alternative models converge to same upper and lower steady states (X_{T }= 95 for upper steady state and X_{T }= 0.0001 for lower steady state Figure 3a. The full model is represented by solid red line (up state) and solid blue line (down state). The various dashed lines present the lower dimensional approximations. We see that even though in condition I the steady state behavior is well approximated by the reduced systems, the dynamics are not. By setting derivatives to zero and replacing some dynamical variables by functions we have produced reduced dynamics that are faster than the full model dynamics. These types of approximations can produce reasonable dynamics if there is a clear separation of time scales, in which case setting the derivatives of the fast dynamics to zero causes only minor differences in the dynamics. Here, although degradation is the slowest dynamical variable, other processes such as dephosphorylation can also be slow, and therefore there is no clear time scale separation. For condition II (Figure 3b) even the steady state behavior of the systems is not well approximated by the reduced equations.
Next, we analyze the dependence of the steady state solution of this loop on different system parameters. This is similar to what we do in Figure 2, but here we show the full bifurcation diagrams resulting from a change in parameters. We also compare cases where the approximations used are appropriate or not. First, we select the degradation rate (λ) of αCaMKII as a bifurcation parameter. Here, degradation represents the removal of αCaMKII either by general cellular degradation pathway or by diffusion from a certain specific location e.g., active synapses. We analyzed the effect of degradation parameter on the characteristics of polyadenylation loop under the two different set of parameters (conditions I and II). We show results of this analysis in Figure 4. On the left we show the G' and F functions (Figure 4a,c), and on the right the complete bifurcation diagrams in terms of λ (Figure 4b,d). Figures 4 a, b are for condition I and 4 c, d is for condition II.
Figure 4. Steady state solution characteristics of molecular loop with respect to parameter λ. This contain an analytical solution (a, c) and numerical/analytical bifurcation diagrams (b, d). Approximate analytical solution is developed through graphing equation 6, whereas the numerical/analytical bifurcation diagrams are developed through tracking the steady state behavior of all three models with respect to the degradation rate. Two different set of parameters are compared. For first set of parameters the approximate solution and numerical/analytical bifurcation diagrams from three alternative models converge to same upper, stable steady state solution branch (a, b). For the second set of parameters the approximate solution and numerical/analytical bifurcation diagrams from three alternative models does not converge to same upper, stable steady state and solution branch(c, d). The approximate analytical solution for first set of parameters (a) is located at two different degradation rates (The solid blue line#"1" represent λ = 0.0001s^{1 }and other solid blue line#"2" represent λ = 0.0003 s^{1 }). This method locates the upper stable steady state at X_{T }= 95, while the lower and unstable steady state are at X_{T }= 0 and X_{T }= 9.4 when degradation rate is set at 0.0001 s^{1 }(curve 1). As degradation rate is increased to 0.0003 s^{1 }(curve 2) the new solution is located at X_{T }= 31 (upper stable steady state), X_{T }= 0.0001 (lower stable steady state) and X_{T }= 9.2 (unstable steady state). The numerical and analytical bifurcation diagrams for first set of parameters (b) are developed through tracking the steady state behavior of all three models. Three bifurcation diagrams exactly match with each other for the entire range of degradation parameter (dotted blue line represent the bifurcation diagram based on full model, whereas solid blue lines represent the bifurcation diagrams based on reduced three differential equation and a single equation model). The approximate analytical solution for second set of parameters (c) is also located at two different degradation rates (The solid blue curve "1" represent λ = 0.0001s^{1}, and other solid blue curve "2" represent λ = 0.0003 s^{1}). This method locates the upper stable steady state at X_{T }= 95, while the lower and unstable steady state are at X_{T }= 0.0001 and X_{T }= 9.4 when degradation rate is set at 0.0001s^{1 }(curve 1). As degradation rate is increased to 0.0003 s^{1}(curve 2) the new solution is located at X_{T }= 31 (upper stable steady state), X_{T }= 0.0001 (lower stable steady state) and X_{T }= 9.2 (unstable steady state). The numerical and analytical bifurcation diagrams for second set of parameters (d) are developed. The bifurcation diagram from full model does not match with the bifurcation diagrams from reduced and single equation model.
We show the F functions for two different values of λ, where "1" denotes λ = 0.0001s^{1 }and "2" denotes λ = 0.0003s^{1}. By this method we can graphically locate the values of X^{T }at the fixed points. For condition I the upper stable steady state at X_{T }= 95, while the lower and unstable steady state are at X_{T }= 0.0001 and X_{T }= 9.4 when degradation rate is set at 0.0001 s^{1 }(curve "1", Figure 4a). As the degradation rate is increased to 0.0003 s^{1 }(curve "2", Figure 4a) the new solution is located at X_{T }= 31 (upper stable steady state), X_{T }= 0 (lower stable steady state) and X_{T }= 9.2 (unstable steady state). Numerical bifurcation diagrams (Figure 4b) are developed through tracking the steady state behavior of all three levels of simplification with respect to the degradation rate of αCaMKII. The solutions contained in these bifurcation diagrams are tracked through a numerical bifurcation package [30], for the 6 and three dimensional models, by simply finding the crossover in the one dimensional model, and by finding zeros of a polynomial in the polynomial approximation. The four bifurcation diagrams are nearly identical for the entire range of degradation parameter. In contrast, to these results when simulations are carried out with second set of parameters (condition II) the steady state solution characteristics of full scale model (Figure 4ddotted blue line represent the bifurcation diagram based on full model) does not match either with the three state reduced model or single equation model or the polynomial model (Figure 4dsolid blue lines represent the bifurcation diagrams based on reduced three differential equation and a single equation model). What these results also indicate is that the most significant approximation made here is in equation 1, and it fails when the conditions of this approximation are not met.
We also analyzed the steady state solution characteristics of this molecular switch with respect to CPEB1 activation parameter k_{5}. First, the simulations are implemented with parameters from condition I (Figure 5, a, b) and then with condition II (Figure 5, c, d). Here, again the steady state solution characteristics of the full scale model do not match (Figure 5ddotted blue line represent the bifurcation diagram based on full model, whereas solid blue lines represent the bifurcation diagrams based on reduced three differential equation and a single equation model) with reduced versions when condition two is implemented, however, with parameters representing the condition one there is a complete match (Figure 5bdotted blue line represent the bifurcation diagram based on full model, whereas solid blue lines represent the bifurcation diagrams based on reduced three differential equation and a single equation model).
Figure 5. Steady state solution characteristics of molecular loop with respect to CPEB1 activation parameter k_{5}. The parametric steady state solution characteristics are developed through an approximate analytical solution (a, c) and numerical/analytical bifurcation diagrams (b, d). The approximate analytical solution is developed through graphing equation 6, whereas the numerical/analytical bifurcation diagrams are developed through tracking the steady state behavior of all three models with respect to k_{5}. Two different set of parameters are compared. For first set of parameters the approximate solution and numerical/analytical bifurcation diagrams from three alternative models converge to same upper, stable steady state solution branch (a, b). For the second set of parameters the approximate solution and numerical/analytical bifurcation diagrams from three alternative models does not converge to same upper, stable steady state and solution branch (c, d). The approximate analytical solution for first set of parameters (a) is located at two different values of rate constant k5 (The dotted red curve "1" represent k5 = 0.0072 μM^{1}.s^{1}, second dotted red curve "2" represent k5 = 0.0001 μM^{1 }.s^{1 }). This method locates the upper stable steady state at X_{T }= 95, while the lower and unstable steady state are at X_{T }= 0.0001 and X_{T }= 9.4 when k_{5 }is set at 0.0072μM^{1 }.s^{1 }(curve "1"). As k_{5 }is decreased to 0.0001 new solution is located at X_{T }= 70 (upper stable steady state), X_{T }= 0.0001 (lower stable steady state) and X_{T }= 16 (unstable steady state). The numerical and analytical bifurcation diagrams for first set of parameters (b) are developed through tracking the steady state behavior of all three models with respect to k5. The three bifurcation diagrams are exactly matching with each other for the entire range of activation parameter (dotted blue line represent the bifurcation diagram based on full model, whereas solid blue lines represent the bifurcation diagrams based on reduced three differential equation and a single equation model). The approximate analytical solution for second set of parameters (c) is also located at two different values of rate constant k_{5 }(The dotted red curve "1" represent k_{5 }= 0.0072 μM^{1}.s^{1}, second dotted red curve "2" represent k5 = 0.0001 μM^{1}.s^{1}). The numerical and analytical bifurcation diagrams for second set of parameters (d) are developed through tracking the steady state behavior of all three models with respect to k5. Full model bifurcation diagram does not match with reduced and single equation model.
Discussion
We have postulated that a translationactivation loop can form a bistable switch that could explain the mechanism for maintaining long term memories. Here, we explore a specific case of this hypothesis composed of a positive feedback loop between αCaMKII and its translation regulator CPEB1. The possibility of a positive feedback loop is confirmed by two previous experimental observations (a) phosphorylation of CPEB1 regulates the synthesis of αCaMKII molecules through polyadenylation of αCaMKIImRNA [20] (b) αCaMKII phosphorylates the CPEB1 molecule [25,26]. Thus, both the αCaMKII and CPEB1 interact through a closed loop. Based on the elementary enzyme kinetics (Method A) the dynamics of this loop can be characterized through a high dimensional system of ordinary differential equations. One can study the characteristics of such a system by numerically integrating these differential equations and with different initial conditions and kinetic rate parameters, and extensively search for the possibility of a two distinct stable steady states. This procedure itself is very tedious, time consuming and inaccurate.
Here, we develop a systematic approach to reduce the dimensionality of a high dimensional ODE based model of αCaMKIICPEB1 molecular loop. Our simplification strategy is devised such that the essential features of full scale model are preserved. This simplification process is based on introducing the two key approximations. 1. The synthesisdegradation curve (equation 1) and 2. The fitting function curve (equation 7). Through these approximations the full scale model is reduced to three lower dimensional models. (a) A three dimensional model comprising of equations 6, B6, and B12. (b) A one dimensional model (equation 6), (c) A one dimensional fourth order polynomial based model (B18). We can extract the fixed point solution characteristics of this system from all three reduced models however, only reduced model (a) and (b) can provide the dynamical characteristics of this system. Since, during the simplification process it is possible to loose some features therefore, it would be logical to compare the performance of these three reduced models with a full scale model of this loop.
Our results show that the key to performance of three reduced models in locating fixed points of the system is the approximation introduced through equation 1. Since, in this approximation the CaMKII bound in biochemical complexes also degrade along with free CaMKII, which is in contrast to full scale model. Thus, for parametric conditions where the amount of bound CaMKII is negligible the three reduced models yield matching results to full scale model. However, for conditions where there is a nonnegligible amount of bound CaMKII the results of three models are not matching with full model. Interestingly, however the results of these three reduced models are matching with each other. Here, this phenomenon of CaMKII trapping in biochemical complexes is simply modeled through an autophosphorylation parameter k_{3}. For the k_{3 }set at lower values the amount of CaMKII bound is large and thus leading to disagreements between the results from full scale model and reduced models.
The analysis of this paper also provides some intuition on the complex process of polyadenylation based translation. Here, through a systematic simplification process we develop a single variable dynamical equation (eq. 6) which shows that the change in total CaMKII concentration is due to a balance between protein synthesis and degradation. At steady state this equation is composed of two terms i.e., synthesis and degradation term. Thus, the polyadenylation based switching is essentially due to a balance between the synthesis of new proteins and degradation of old ones. If the amount of synthesis is at much faster rate then the degradation rate the balance will shift towards single state high protein concentration, similarly if the degradation rate is much faster then new protein synthesis the balance will tip towards single lower state of switch. So for reversible switching a balance should be maintained between protein synthesis and degradation. This equation also provides an easy method to analyze the effect of other system variables on switching characteristics. For example one parameter which could critically affect the polyadenylation switching is the CPEB1 activation rate through CaMKII i.e., k_{5}. For larger value of this parameter the new synthesis rate will increase and at lower value the new synthesis rate will decrease.
In this paper, we simplify the high dimensional ODE model of αCaMKIICPEB1 loop in such a way that these simplifications preserve the characteristics properties of the full scale ODE model. By reducing the high dimensional model to the minimal plausible scenario, we were able to obtain a single algebraic expression that characterizes the fixed points of these dynamics. We then used these algebraic equations to develop the analytical bifurcation diagrams by perturbing the various parameters (λ and k5). Using this analytical expression we explain how the loss of αCaMKII in synapses due to protein degradation is balanced by synthesis of new αCaMKII molecules. This analysis explains how the pair of two molecules in active and inactive forms and a synthesis based positive feedback loop can lead to a bistable switch. On the basis of these bifurcation diagrams we show that even the simplified version of a polyadenylation based model of αCaMKII can exhibits an up and down state and by perturbing these parameters the system can toggle between an up and a down state. The activity induced increase in total amount of αCaMKII can be maintained for long period of time due to a balance between the new synthesis and degradation of αCaMKII molecule. As we set equation 1 equals to zero the resulting equation 2 describe the steady state value of total amount of αCaMKII. It also shows that total concentration of CaMKII is directly proportional to the amount of CPEB1 phosphorylated and is inversely proportional to the degradation rate. Therefore, when the fraction of CPEB1 phosphorylated is high the amount of CaMKII generated through polyadenylation increases, which will out balance the amount of αCaMKII removed through degradation and there is a net upregulation of αCaMKII concentration. Similarly this equation also explains that if degradation rate is too fast the compensation provided by new synthesis of αCaMKII molecules will not be enough to balance the loss of CaMKII due to degradation. Our simplification captures the essential properties of protein translation through polyadenylation. It shows that how the activity induced signal converts the inactive form of αCaMKII into active, which is further amplified through an autophosphorylation loop. The active and phosphorylated αCaMKII in turn drives the synthesis of a new inactive αCaMKII molecule, through phosphorylating the CPEB1. Here, all these steps are analytically proven as shown by the equations (B1B9).
Apart, from the approximation introduced through equation 1 another approximation is introduced into this analysis when we tried to develop a single algebraic equation 6 based version of this system. In equation 6 we developed an approximate function XP = g ( XT). In order, to develop a completely accurate solution this function should be analytically developed, however, for this system it was not possible to invert the highly nonlinear system of XT as a function of XP. Here, we developed the approximate function XP = g (XT) through a fitting function routine in matlab. The approximate fitting function proved to be very sufficient good in approximating both the upper, lower stable steady state solution and unstable steady state solution., however, it could not accurately approximate the lower stable steady state solution.
The numerical values of different parameters of this molecular loop are described in table 1. Many of these parameters are extracted from previous observations based on experimental and simulation work [16,19,20,2534]. For example parameters like αCaMKII degradation rates, calcium/calmodulin binding and unbinding rate, rate parameters for αCaMKII autophosphorylation loop and rate parameters describing the new synthesis of αCaMKII protein are taken from previous experimental and simulation based observations [27,28,3133], and [34]. Some other rate parameters such as activation of CPEB1 through phosphorylated αCaMKII are obtained through parameters scaling and matching the observed experimental dynamics. The parameters of table 1 and bifurcation diagrams of Figure 4 and 5 indicate that this molecular system is very robust. For example the k_{5 }parameter (describing αCaMKII mediated activation of CPEB1) exhibits the bistable character over a range of three orders of magnitude. This means that even a small activation of CPEB1 through αCaMKII will induce the bistable character to this loop and large strength of αCaMKII activity inhibition is required to reverse the upstate to down. This leads to interesting observation where many of αCaMKII inhibitors are not able to reverse the established LLTP when applied during the maintenance phase of memory formation [33]. The reduced model of this work and fullscale model [30] can make these predictions.
Conclusions
This study provides a systematic method to simplify, approximate and analyze a molecular model of polyadenylation loop. This model of polyadenylation loop is based upon denovo synthesis of new proteins & their activation. The polyadenylation loop operates through a positive feedback loop between a protein and its translation factor. This work shows how to extract low dimensional systems that can be used to obtain analytical solutions for the fixed points of the system and to describe the dynamics of the system. The methods used here have general applicability to the formulation and analysis of many molecular networks.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
NA carried out the modeling, simulation and analysis of this work. HS contributed in conceptual design and polishing of this work. All authors read and approved the final manuscript.
References

Jacob F, Monod J: Genetic regulatory mechanisms in the synthesis of proteins.
J Mol Biol 1961, 3:318356. PubMed Abstract  Publisher Full Text

Ozbudak E, Thattai M, Lim H, Shraiman B, van Oudenaarden A: Multistability in the lactose utilization network of Escherichia coli.
Nature 2004, 427:737740. PubMed Abstract  Publisher Full Text

Ferrell JE: Selfperpetuating states in signal transduction: positive feedback, double negative feedback and bistability.
Curr Opin Cell Biol 2002, 14:140148. PubMed Abstract  Publisher Full Text

Angeli D, Ferrell JE, Sontag ED: Detection of multi stability, bifurcations, and hysteresis in a large class of biological positivefeedback systems.
PNAS 2004, 101:7.
18221827

Ferrell JE, Xiong W: Bistability in cell signaling: how to make continuous processes discontinuous, and reversible processes irreversible.
Chaos 2001, 11:227236. PubMed Abstract  Publisher Full Text

Novik A, Wiener M: Enzyme induction as an allornone phenomenon.

Chon M, Horibata K: Inhibition by glucose of the induced synthesis of the βGALACTOSIDEenzyme system of Escerichia coli.

Ptashne M: A Genetic Switch: Phage and Higher Organisms. Blackwell, Oxford; 1992.

Ferrell JE Jr, Machleder EM: The biochemical basis of an allornone cell fate switch in Xenopus Oocytes.

Bagowski CP, Ferrell JE: The JNK cascade as a biochemical switch in mammalian cells: ultrasensitive and allornone responses.
Curr Biol 2001, 11:11761182. PubMed Abstract  Publisher Full Text

Bhalla US, Ram PT, Iyengar R: MAP kinase phosphatase as a locus of flexibility in a mitogenactivated protein kinase signaling network.
Science 2002, 297:10181023. PubMed Abstract  Publisher Full Text

Cross FR, Archambault V, Miller M, Klovstad M: Testing a mathematical model of the yeast cell cycle.
Mol Bio Cell 2002, 13:5270. Publisher Full Text

Sha W, Moore J, Chen K, Lassaletta AD, Yi CS, Tyson JJ, Sible JC: Hysteresis drives cellcycle transitions in xenopus laevis egg extract PNAS.
USA 2003, 100:975980. Publisher Full Text

Pomerening JR, Sontag ED, Ferrel JE Jr: Building a cell cycle oscillator: hysteresis and bistability in the activation of Cdc2.
Nat Cell Biol 2003, 5:346351. PubMed Abstract  Publisher Full Text

Bagowski CP, Besser J, Frey CR, Ferrel JE Jr: The JNK cascade as a biochemical switch in mammalian cells: ultrasensitive and allornone responses. Curr.

Lisman JE: A mechanism for memory storage insensitive to molecular turnover a bistable autophosphorylating kinase.
PNAS 1985, 82:30553057. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Xiong W, Ferrel JE Jr: A positivefeedbackbased 'memory module' that governs a cell fate decision.
Nature 2003, 426:460465. PubMed Abstract  Publisher Full Text

Reynolds AR, Tischer C, Verveer PJ, Rooks O, Bastiaens PI: EGFR activation coupled to inhibition of tyrosine phosphatases causes lateral signal propagation.
Nat Cell Biol 2003, 5:447453. PubMed Abstract  Publisher Full Text

Wells DG, Dong X, Quinlan EM, Huang YS, Bear MF, Richter JD, Fallon JR: A role for cytoplasmic polyadenylation element in NMDA receptorregulated mRNA translation in neurons.
J Neurosci 2001, 21:95419548. PubMed Abstract  Publisher Full Text

Wells DG, Richter JD, Fallon JR: Molecular mechanisms for activityregulated protein synthesis in the synaptodendritic compartment.
Curr Opin Neurobiol 2000, 10:132137. PubMed Abstract  Publisher Full Text

Wilutz CJ, Wormington M, Peltz SW: The CapToTail guide to mRNA turnover.

Gallie DR: The Cap and poly(A) tail function synergistically regulate mRNA translational efficiency.
Genes & Development 1991, 5:21082116. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

De Moor C, Meijer H, Lissenden S: Mechanisms of translational control by the 3'UTR in development and differentiation.
Cell and Development Biology 2005, 16:4958. Publisher Full Text

Piccioni F, Zappavigna V, Verrotti AC: Translational regulation during oogenesis and early development: The cappoly(A) tail relationship. C.R.
Biologies 2005, 328:863881. PubMed Abstract  Publisher Full Text

Atkins CM, Davare MA, Oh M, Derkach V, Soderling TR: Bidircetional regulation of cytoplasmic polyadenylation elementbinding protein phosphorylation by Ca^{+2}/Calmodulindependent kinase II and protein phosphatase1 during Hippocampal Longterm potentiation.
J Neurosci 2005, 25:56045610. PubMed Abstract  Publisher Full Text

Atkins CM, Nozaki N, Shigeri Y, Soderling TR: Cytoplasmic polyadenylation element binding proteindependent protein synthesis is regulated by Calcium/Calmodulindependent protein kinase II.
J Neurosci 2004, 24:51935201. PubMed Abstract  Publisher Full Text

Ehlers MD: Activity level controls postsynaptic composition and signaling via the ubiquitinproteasome system.
Nat Neurosci 2003, 6:231242. PubMed Abstract  Publisher Full Text

Aslam N, Kubota Y, Shouval HZ: Translational switch for longterm maintenance of synaptic plasticity.

Sommese AJ, Wampler WC: The numerical solution of systems of polynomials: Arising in Engineering and Science. Volume 2005. World Scientific Publishing Co N.J;

Dhooge A, Govaerts W, Kuznetsov , Yu A: MATCONT: A MATLAB package for numerical bifurcation analysis of ODE's ACM TOMS.
2003, 29:141164.

Bradshaw JM, Kubota Y, Meyer T, Schulman H: An ultrasensitive Ca+2/calmodulindependent protein kinase IIprotein phosphatase 1 switch facilitates specificity in postsynaptic calcium signaling.
Proc Natl Acad Sci USA 2003, 100:1051210517. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kubota Y, Bower JM: Transient versus asymptotic dynamics of CaMKII: possible role of phosphatase.
J Comput Neurosci 2001, 11:263279. PubMed Abstract  Publisher Full Text

Otmakhov N, Griffith LC, Lisman JE: Postsynaptic inhibitors of calcium/calmodulindependent protein kinase type II block induction but not maintenance of pairinginduced longterm potentiation.
J Neurosci 1997, 17:53575365. PubMed Abstract  Publisher Full Text

Sanhueza M, McIntyre CC, Lisman JE: Reversal of synaptic memory by Ca+2/calmodulindependent protein kinase II inhibitor.
J Neurosci 2007, 27:51905199. PubMed Abstract  Publisher Full Text