Search tips
Search criteria 


Logo of bmcsysbioBioMed Centralsearchsubmit a manuscriptregisterthis articleBMC Systems Biology
BMC Syst Biol. 2012; 6: 12.
Published online 2012 February 15. doi:  10.1186/1752-0509-6-12
PMCID: PMC3359155

Regulation of cytoplasmic polyadenylation can generate a bistable switch



Translation efficiency of certain mRNAs can be regulated through a cytoplasmic polyadenylation process at the pre-initiation 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/Calmodulin-dependent protein kinase II (αCaMKII) translation through the phosphorylation/dephosphorylation cycle of CPEB1.


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 non-linear ordinary differential equations can describe the dynamic characteristics of the polyadenylation loop. Here, we simplified this high-dimensional 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.


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.


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 non-linear interactions to exhibit a bistable behavior. Previous, experimental work has described several examples of naturally occurring bistable system [6-18]. 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., activation-deactivation 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 αCaMKII-mRNA 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 self-sustained 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 αCaMKII-CPEB1 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.


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 Michalis-Menten 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 pseudo-steady state assumptions the differential equations representing the complexes [C1-C3] can be eliminated:


The unphosphorylated CPEB1 (Y) is related to phosphorylated CPEB1 (YP).


Where, YT 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 (YP) from R1-R6


By requiring a steady state dYPdt=0 and substituting a=k6*k5(k55+k6) we get


Also from (B-1,B-2, B-3, B-4, B-5 and B-7) we can obtain the value of X.


Where, N and M are constants defined as N=k5(k55+k6),M=k1(k2+k3)

As described in the result section the equation 2 is as follows:


Substituting C3 from (B-3) in [2]


Placing YP from (B-7) into (B-9)


Where, P7 = P*k7 and P4 = P*k4. Defining, U = (Ca+2)4-CaM, b=k9*k8*T(k88+k9) and c = b*YT or c=YT*k9*k8*T(k88+k9)

We obtain


The equation 1 describes the rate of change of total concentration of αCaMKII (XT) 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 R1-R6 an equation representing the dynamics of phosphrylated αCaMKII can be constructed.


The equation B-12 is further simplified in only two variables i.e., XT and XP by placing the dYPdt=0 and putting the values of X, C1, C2, Y, and YP from B-1 to B-8 and further simplifying we get the following expression.


The I1,I2,I3 and I4 are defined as follows:


The expression in B-14 to B-17 are defined as k7 = k7 *P, k4 = k4 *P,a1 = (k3 *k1 )./(k2 + k3 ) a2 = k10 *U, DI = k4 + k1010 , d1 = (k6 .*k5)./(k55 + k6)

The expression of B-13 gives us a function where XT = f(XP), It is almost impossible to invert this function to get XP = f-1(XT)so we numerically approximated it through function fitting.

Equation B-12 with the aid of equation B1-B 13 can be transformed to following polynomial in terms of XP.


Where, the coefficients z12,z13 and z14 are defined by following expressions (B-19-B-32). The equation (B-18) is obtained through further simplifying (B-13) such that a fourth order polynomial is obtained in terms of XP. This 4'th order polynomial has an analytical solution because it has no zero order term and therefore has one solution XP = 0, and the other solutions are the solutions of a third order polynomial.


Analysis of the model of αCaMKII synthesis through polyadenylation

Our simplified model of a self-sustained polyadenylation of αCaMKII-mRNA (Figure (Figure1)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 αCaMKII-mRNA (Figure (Figure1)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 pre-initiation 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 self-sustaining αCaMKII-CPEB 1 molecular loop. The αCaMKII, molecule can be inactive, or active and phosphorylated state, while CPEB1 can be in unphosphorylated ...

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 XP.

Similarly CPEB is denoted as Y and phosphorylated CPEB1 as YP. 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 two-state assumption of the αCaMKII is a significant simplification compared to more realistic multi-state 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 (C1-C3) to zero (Method B), to obtain a pseudo steady state approximation, similar to that used in the standard michaelis-menten approximation. This reduces our system to 3 dynamical equations, and three functional expressions, as described by equations 1, B-6, and B-12. 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, C3, as described by equation B-3 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 k9 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 XT 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.

Defining dXTdt=H(C3,XT) we have that:


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 YT, 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 Km1/2 of this process is P7/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 Xp to XT, i.e., XP = g (XT), 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 XT as a function of XP [(Method B equation 13)], it is difficult to analytically invert it to obtain an analytical expression for XP = g ( XT). 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 S1.The solid black line in additional file 1, Figure S1 represents the original function {Method B B-13}, 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:


Here, the XP is obtained from fitting function 7 and F (XT) 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 (Figure2a).2a). The function G is shown as dotted red line and F as solid blue line. This graphical representation (Figure (Figure2a)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 XT there is hardly any synthesis, and degradation dominates, so the system converges to the low fixed point at zero. At higher values of XT 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 XT increases until it reaches the third intersection, which is the upper stable fixed point. Due to the saturation of G, at higher levels of XT, 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 αCaMKII-CPEB1 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 ...

Next we show the impact of changing some of the system parameters on the steady state solution characteristics of equation 6 (Figure (Figure2b2b and and2c).2c). Changing the degradation rate only affects the degradation curve (Figure (Figure2bsolid2bsolid blue line). As the degradation rate is increased from slower (Figure (Figure2bsolid2bsolid blue line 1) to much faster (Figure (Figure2bsolid2bsolid blue line 4) the system moves from robustly bistable to a mono-stable 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 (Figure2bsolid2bsolid 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 (Figure2bsolid2bsolid blue lines 2-4). As the degradation rate is decreased, and the system as shifts from mono-stable to bistable, the lower fixed point does not change, the value of XT at the unstable fixed point is increased, and the value of XT at the upper fixed point first rapidly increases and then plateaus (Figure (Figure2bsolid2bsolid blue line 2,3,4).

Other parameters affect only the synthesis curve (G'). One such example is parameter k5, which quantifies the CaMKII mediated activation of CPEB1 (Figure (Figure2c).2c). When this activation parameter is set at low value (Figure (Figure2cdotted2cdotted 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 (Figure2cdotted2cdotted red line 3,2,1) and shifting the system from mono-stable 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 z12-z14 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 synthesis-degradation loop (Figure (Figure1)1) and the various approximations used to simplify it. We can integrate the dynamics either using the full system of ODE's (Method A, A1-A6) or by using the reduced three dimensional systems (Method B, equations 6, B-6, 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, A1-A6), the reduced three dimensional system (6, B-6, 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 (Figure3)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 (Figure3a),3a), whereas the second set of parameters (Condition II) represents a case where this approximation does not hold (Figure (Figure3b).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 auto-phosphorylation parameter k3. For condition I we set k3 = 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 (Table1).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 up-state.

Table 1
The numerical values of different parameters of bistable molecular loop.
Figure 3
The dynamic characteristics of aCaMKII-CPEB 1 molecular loop. These characteristics are developed through three alternative models. First model is based on differential equations A1-A6 (Full model, represented by solid red line (up state) and solid blue ...

For condition I the dynamics from all three alternative models converge to same upper and lower steady states (XT = 95 for upper steady state and XT = 0.0001 for lower steady state Figure Figure3a.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 (Figure3b)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 Figure2,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 Figure4.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 ...

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 XT at the fixed points. For condition I the upper stable steady state at XT = 95, while the lower and unstable steady state are at XT = 0.0001 and XT = 9.4 when degradation rate is set at 0.0001 s-1 (curve "1", Figure Figure4a).4a). As the degradation rate is increased to 0.0003 s-1 (curve "2", Figure Figure4a)4a) the new solution is located at XT = 31 (upper stable steady state), XT = 0 (lower stable steady state) and XT = 9.2 (unstable steady state). Numerical bifurcation diagrams (Figure (Figure4b)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 cross-over 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 (Figure4ddotted4ddotted 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 (Figure4dsolid4dsolid 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 k5. 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 (Figure5ddotted5ddotted 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 (Figure5bdotted5bdotted 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 k5. The parametric steady state solution characteristics are developed through an approximate analytical solution (a, c) and numerical/analytical bifurcation ...


We have postulated that a translation-activation 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 αCaMKII-mRNA [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 αCaMKII-CPEB1 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 synthesis-degradation 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, B-6, and B12. (b) A one dimensional model (equation 6), (c) A one dimensional fourth order polynomial based model (B-18). 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 non-negligible 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 k3. For the k3 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., k5. 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 αCaMKII-CPEB1 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 up-regulation 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 auto-phosphorylation 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 (B1-B9).

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 non-linear 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 table1.1. Many of these parameters are extracted from previous observations based on experimental and simulation work [16,19,20,25-34]. For example parameters like αCaMKII degradation rates, calcium/calmodulin binding and un-binding rate, rate parameters for αCaMKII auto-phosphorylation loop and rate parameters describing the new synthesis of αCaMKII protein are taken from previous experimental and simulation based observations [27,28,31-33], 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 table11 and bifurcation diagrams of Figure Figure44 and and55 indicate that this molecular system is very robust. For example the k5 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 up-state to down. This leads to interesting observation where many of αCaMKII inhibitors are not able to reverse the established L-LTP when applied during the maintenance phase of memory formation [33]. The reduced model of this work and full-scale model [30] can make these predictions.


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 de-novo 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.

Supplementary Material

Additional file 1:

Figure S1 The fitting performance of equation7 with original function B-13 (Method B).


  • Jacob F, Monod J. Genetic regulatory mechanisms in the synthesis of proteins. J Mol Biol. 1961;3:318–356. doi: 10.1016/S0022-2836(61)80072-7. [PubMed] [Cross Ref]
  • Ozbudak E, Thattai M, Lim H, Shraiman B, van Oudenaarden A. Multistability in the lactose utilization network of Escherichia coli. Nature. 2004;427:737–740. doi: 10.1038/nature02298. [PubMed] [Cross Ref]
  • Ferrell JE. Self-perpetuating states in signal transduction: positive feedback, double negative feedback and bistability. Curr Opin Cell Biol. 2002;14:140–148. doi: 10.1016/S0955-0674(02)00314-9. [PubMed] [Cross Ref]
  • Angeli D, Ferrell JE, Sontag ED. Detection of multi stability, bifurcations, and hysteresis in a large class of biological positive-feedback systems. PNAS. 2004;101:7. 1822-1827. [PubMed]
  • Ferrell JE, Xiong W. Bistability in cell signaling: how to make continuous processes discontinuous, and reversible processes irreversible. Chaos. 2001;11:227–236. doi: 10.1063/1.1349894. [PubMed] [Cross Ref]
  • Novik A, Wiener M. Enzyme induction as an all-or-none phenomenon. PNAS USA. 1975;43:553–566. [PubMed]
  • Chon M, Horibata K. Inhibition by glucose of the induced synthesis of the β-GALACTOSIDE-enzyme system of Escerichia coli. Analysis of Maintenance J Bacteriol. 1959;78:601–612. [PMC free article] [PubMed]
  • Ptashne M. A Genetic Switch: Phage and Higher Organisms. Blackwell, Oxford; 1992.
  • Ferrell JE Jr, Machleder EM. The biochemical basis of an all-or-none cell fate switch in Xenopus Oocytes. Science. 1998;297:1018–1023. [PubMed]
  • Bagowski CP, Ferrell JE. The JNK cascade as a biochemical switch in mammalian cells: ultrasensitive and all-or-none responses. Curr Biol. 2001;11:1176–1182. doi: 10.1016/S0960-9822(01)00330-X. [PubMed] [Cross Ref]
  • Bhalla US, Ram PT, Iyengar R. MAP kinase phosphatase as a locus of flexibility in a mitogen-activated protein kinase signaling network. Science. 2002;297:1018–1023. doi: 10.1126/science.1068873. [PubMed] [Cross Ref]
  • Cross FR, Archambault V, Miller M, Klovstad M. Testing a mathematical model of the yeast cell cycle. Mol Bio Cell. 2002;13:52–70. doi: 10.1091/mbc.01-05-0265. [PMC free article] [PubMed] [Cross Ref]
  • Sha W, Moore J, Chen K, Lassaletta AD, Yi CS, Tyson JJ, Sible JC. Hysteresis drives cell-cycle transitions in xenopus laevis egg extract PNAS. USA. 2003;100:975–980. doi: 10.1073/pnas.0235349100. [PubMed] [Cross Ref]
  • 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:346–351. doi: 10.1038/ncb954. [PubMed] [Cross Ref]
  • Bagowski CP, Besser J, Frey CR, Ferrel JE Jr. The JNK cascade as a biochemical switch in mammalian cells: ultrasensitive and all-or-none responses. Curr. Biol. 2003;13:315–320. [PubMed]
  • Lisman JE. A mechanism for memory storage insensitive to molecular turnover a bistable autophosphorylating kinase. PNAS. 1985;82:3055–3057. doi: 10.1073/pnas.82.9.3055. [PubMed] [Cross Ref]
  • Xiong W, Ferrel JE Jr. A positive-feedback-based 'memory module' that governs a cell fate decision. Nature. 2003;426:460–465. doi: 10.1038/nature02089. [PubMed] [Cross Ref]
  • 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:447–453. doi: 10.1038/ncb981. [PubMed] [Cross Ref]
  • Wells DG, Dong X, Quinlan EM, Huang YS, Bear MF, Richter JD, Fallon JR. A role for cytoplasmic polyadenylation element in NMDA receptor-regulated mRNA translation in neurons. J Neurosci. 2001;21:9541–9548. [PubMed]
  • Wells DG, Richter JD, Fallon JR. Molecular mechanisms for activity-regulated protein synthesis in the synapto-dendritic compartment. Curr Opin Neurobiol. 2000;10:132–137. doi: 10.1016/S0959-4388(99)00050-1. [PubMed] [Cross Ref]
  • Wilutz CJ, Wormington M, Peltz SW. The Cap-To-Tail guide to mRNA turnover. Nature. 2001;2:237–246. [PubMed]
  • Gallie DR. The Cap and poly(A) tail function synergistically regulate mRNA translational efficiency. Genes & Development. 1991;5:2108–2116. doi: 10.1101/gad.5.11.2108. [PubMed] [Cross Ref]
  • 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:49–58. doi: 10.1016/j.semcdb.2004.11.007. [PubMed] [Cross Ref]
  • Piccioni F, Zappavigna V, Verrotti AC. Translational regulation during oogenesis and early development: The cap-poly(A) tail relationship. C.R. Biologies. 2005;328:863–881. doi: 10.1016/j.crvi.2005.05.006. [PubMed] [Cross Ref]
  • Atkins CM, Davare MA, Oh M, Derkach V, Soderling TR. Bidircetional regulation of cytoplasmic polyadenylation element-binding protein phosphorylation by Ca+2/Calmodulin-dependent kinase II and protein phosphatase-1 during Hippocampal Long-term potentiation. J Neurosci. 2005;25:5604–5610. doi: 10.1523/JNEUROSCI.5051-04.2005. [PubMed] [Cross Ref]
  • Atkins CM, Nozaki N, Shigeri Y, Soderling TR. Cytoplasmic polyadenylation element binding protein-dependent protein synthesis is regulated by Calcium/Calmodulin-dependent protein kinase II. J Neurosci. 2004;24:5193–5201. doi: 10.1523/JNEUROSCI.0854-04.2004. [PubMed] [Cross Ref]
  • Ehlers MD. Activity level controls postsynaptic composition and signaling via the ubiquitin-proteasome system. Nat Neurosci. 2003;6:231–242. doi: 10.1038/nn1013. [PubMed] [Cross Ref]
  • Aslam N, Kubota Y, Shouval HZ. Translational switch for long-term maintenance of synaptic plasticity. Molecular Systems Biology. 2009. [PMC free article] [PubMed]
  • Sommese AJ, Wampler WC. The numerical solution of systems of polynomials: Arising in Engineering and Science. Vol. 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. pp. 141–164.
  • Bradshaw JM, Kubota Y, Meyer T, Schulman H. An ultrasensitive Ca+2/calmodulin-dependent protein kinase II-protein phosphatase 1 switch facilitates specificity in postsynaptic calcium signaling. Proc Natl Acad Sci USA. 2003;100:10512–10517. doi: 10.1073/pnas.1932759100. [PubMed] [Cross Ref]
  • Kubota Y, Bower JM. Transient versus asymptotic dynamics of CaMKII: possible role of phosphatase. J Comput Neurosci. 2001;11:263–279. doi: 10.1023/A:1013727331979. [PubMed] [Cross Ref]
  • Otmakhov N, Griffith LC, Lisman JE. Postsynaptic inhibitors of calcium/calmodulin-dependent protein kinase type II block induction but not maintenance of pairing-induced long-term potentiation. J Neurosci. 1997;17:5357–5365. [PubMed]
  • Sanhueza M, McIntyre CC, Lisman JE. Reversal of synaptic memory by Ca+2/calmodulin-dependent protein kinase II inhibitor. J Neurosci. 2007;27:5190–5199. doi: 10.1523/JNEUROSCI.5049-06.2007. [PubMed] [Cross Ref]

Articles from BMC Systems Biology are provided here courtesy of BioMed Central