Analysis of the model of αCaMKII synthesis through polyadenylation
Our simplified model of a self-sustained polyadenylation of αCaMKII-mRNA (Figure ) 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 ) 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.
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 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 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
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 S
1.The solid black line in additional file
1, Figure S
1 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 ). The function G is shown as dotted red line and F as solid blue line. This graphical representation (Figure ) 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.
Next we show the impact of changing some of the system parameters on the steady state solution characteristics of equation 6 (Figure and ). Changing the degradation rate only affects the degradation curve (Figure solid blue line). As the degradation rate is increased from slower (Figure solid blue line 1) to much faster (Figure solid 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 solid 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 solid 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 solid 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 ). When this activation parameter is set at low value (Figure dotted 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 dotted 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 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 synthesis-degradation loop (Figure ) 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 ) 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 ), whereas the second set of parameters (Condition II) represents a case where this approximation does not hold (Figure ). 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 ). 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 1The numerical values of different parameters of bistable molecular loop. |
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 . 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 ) 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 , 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 . On the left we show the G' and F functions (Figure ), and on the right the complete bifurcation diagrams in terms of λ (Figure ). Figures are for condition I and 4 c, d is for condition II.
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 ). As the degradation rate is increased to 0.0003 s
-1 (
curve "2", Figure ) 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 ) 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
dotted 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
solid 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 ) and then with condition II (Figure ). Here, again the steady state solution characteristics of the full scale model do not match (Figure 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) with reduced versions when condition two is implemented, however, with parameters representing the condition one there is a complete match (Figure 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).
1 and Harel Z Shouval1
