MATHEMATICAL ANALYSIS OF A MODULAR NETWORK COORDINATING THE CELL CYCLE AND APOPTOSIS

The cell-division cycle and apoptosis are key cellular processes deregulated during carcinogenesis. Recent work of Aguda and Algar suggests a modular organization of regulatory molecular pathways linking the cellular processes of division and apoptosis. We carry out a detailed mathematical analysis of the Aguda-Algar model to unravel the dynamics implicit in the model structure. In addition, we further explore model parameters that control the bifurcations corresponding to the aforementioned cellular state tran- sitions. We show that this simple model predicts interesting behavior, such as hysteretic oscillations and different conditions in which apoptosis is triggered.


Introduction.
It is now widely accepted among cancer researchers that the cell-division cycle and the cell-death program (called apoptosis) are key cellular processes deregulated during carcinogenesis.Enhanced cell proliferation can be due to increased cell cycling, inhibited apoptosis, or both.Recently, Aguda and Algar [1] reviewed and analyzed the regulatory molecular pathways linking the initiation of the cell cycle and apoptosis.These authors suggested a modular organization of the complex networks and presented a corresponding kinetic model for the cellular state transitions from quiescence (nondividing) to cell cycling, and eventually to apoptosis as the strength of extracellular signaling increases.This kinetic model, however, was not analyzed in detail in the original paper.Here, we extend the mathematical analysis to show that this simple model predicts other interesting behavior, such as hysteretic oscillations and different conditions in which apoptosis is triggered.
The modular structure of the Aguda-Algar model is shown in Figure 1.The lines that end with " " indicate that one module inhibits the other, and the lines that end with "→" indicate that one module activates the other.The modules involved correspond to intracellular signaling, cell cycle, apoptosis, and a control node of transcription factors that stimulate the expression of cell cycle and apoptosis genes.The key transcription factors are suggested to be the proteins E2F-1 and c-Myc.The primary components of the cell-cycle module are the enzymes called cyclindependent kinases, and the key players in the apoptosis module are the caspases.An example of a signaling pathway that impinges on the activation of E2F-1 is the Ras/MAPK pathway which is a cascade of activation involving protein kinases [5].This signaling pathway has been shown in [5] to possess an ultrasensitive switch-like behavior, as depicted in the topmost box of Figure 1.[1] reviewed the literature that points to the proteins E2F-1 and c-Myc as the key transcription factors in this control node.Figure 13 in appendix B gives a summary of the important pathways implicating E2F-1 and c-Myc, and offers details and justification of the modular structure shown in Figure 1 (Fig. 13 is Fig. 3 of [1]).

Aguda and Algar
The interaction labeled "a" in Figure 1 represents the fact that activated cyclindependent kinases (e.g., CDK2) needed for entry into the cell cycle would further enhance the transcriptional activity of E2F-1 via phosphorylation of the retinoblastoma (Rb) protein.The interaction labeled "b" in Figure 1 represents the suggested negative feedback between the Rb/E2F pathway and the Ras signaling pathway ( [6]).The interaction labeled "b'" involves the negative regulation of cyclin A/CDK2 toward E2F-1 through the phosphorylation of the latter's partner, the DP proteins, which are required for binding DNA.Finally, the interaction labeled "c" includes the cell-cycle-enhanced apoptosis via the indirect regulation of caspases by cyclin-dependent kinases (e.g., activation of a caspase by a D-type cyclin, as reported by Mendelsohn et al. [7]).More details and justification of these modulemodule interactions are discussed in [1].
2. Description of the model.A kinetic implementation of the modular model in Figure 1 is shown in Figure 2. The meaning of the notations used in Figure 2 is We will use the same species symbols for their corresponding concentrations.Then, according to either mass-action or Michaelis-Menten kinetics, we get the following system of differential equations: Except for the time-scaling factor ε, equations (1)-( 5) are identical to those used in [1].The cyclic reactions shown in Figure 2-namely, those between the pairs of species (S 1 , S 2 ), (C 1 , C 2 ), and (A 1 , A 2 )-imply that the total concentrations of each pair of species is constant.Of course, the corresponding molecular species are, in general, not constant inside the cell; however, the primary purpose of this simplification is merely to create an ultrasensitive switch (see the two bottom boxes of Fig. 1) that can be modeled using the zeroth-order ultrasensitivity property of cyclic enzyme reactions when the Michaelis constants are close to zero [4].The individual reaction steps in these cyclic reactions are assumed to possess the Michaelis-Menten type of kinetics; that is, the reaction rate has the form vX/(K M + X), where X is the substrate (reactant), v is the maximum reaction rate, and K M is the Michaelis constant.The maximum rate v is itself proportional to the total enzyme concentration catalyzing the step; for example, for the conversion of S 1 to S 2 , v is proportional to Σ, since the latter is the catalyst for this step.In (1) and (3) the interactions take place by mass-action law and degradation.For example, the term 3) is due to degradation promoted by C 2 , and the term k m2 G 2 is due to intrinsic degradation.
We include the various powers of ε to account for the fact that the various processes have different time scales.For negative powers of ε, the smaller the value of ε, the faster the process approaches equilibrium.For positive powers, small ε means that the evolving variable (e.g.A 2 ) changes very slowly.Thus, for example, in the subsystem ( 2)-( 4), the variables S 2 , G 2 , C 2 are fast evolving, and G 2 and S 2 evolve faster than C 2 , as suggested by the parameter values in [1].
We fix We take the parameter reference values (PRVs) to be the following: For these PRVs, the time scale is such that one cell cycle typically takes 50-100 time units.
Except for ε, all other parameters are of the same order of magnitude as those chosen in [1].Our aim is to show that the model is sufficient to account for the following experimental observations (see [2,3,10]): Assertion A. If the signal rate k s is small then the cell is quiescent; if k s is of intermediate size then the cell is cycling; and if k s is large then the cell goes into apoptosis.
We assume that there exists a threshold on the apoptosis marker beyond which apoptosis occurs irreversibly; without loss of generality, we suppose that this threshold is 1  2 .(The parameters were chosen so that the ultrasensitive switch occurs around this value of A 2 ).
Note that, for small ε, the variable G 2 reaches equilibrium much faster than Σ, S 2 , C 2 , and A 2 .Thus we may suppose that, in the time it takes the variable G 2 to reach equilibrium, the change in the other variables is negligible.Therefore, we may equate the right-hand side of (3) to zero and determine G 2 as a function of S 2 and C 2 .Similar fast-slow variable considerations apply to the remaining equations, and they help us determine some rate parameters in our computations as we work with the full dynamical system.
3. Hysteresis loop for signal/cycle interplay.We first consider only the dynamics of S 2 and G 2 which are the fastest variables in the system (1)-( 4) for small ε.Basically, the S 2 -G 2 interaction represents the activation of transcription factors that involves fast protein-protein interactions.The other variables (Σ, C 2 , and A 2 ) are all considered slow, because they are assumed to represent the dynamics of much larger regulatory networks.We equate the right-hand side of (3) to zero, and we solve for G 2 as a function of S 2 and C 2 .Then we equate the right-hand side of (2) to zero and solve for S 2 as a function of Σ.Note that although in general this last equation has two solutions, only the one given by the Goldbeter-Koshland function is biologically relevant (the other solution would lead to a negative value for S 1 , see [4,9]).
Then (4) reduces to the equation More details about the computation of Ψ are given in appendix A. Equations ( 1) and ( 6) capture the dynamics of ( 1)-( 4) restricted to the manifold obtained by making the right-hand sides of ( 2) and (3) equal to zero, and they represent a good approximation of the dynamics of ( 1)-( 4) for small ε.
The S-shaped C 2 -nullcline in Figure 3 is the curve Ψ(C 2 , Σ) = 0; that is, for each fixed Σ, we compute C 2 by solving (2), (3), and (4) at equilibrium.The  upper and lower branches of the S-shaped curve are stable, and the middle branch is unstable.This suggests that if the nullcline of (1) intersects only the middle branch, then as Σ increases to some critical value Σ 2 , there is an upward jump (which can be interpreted as having crossed the checkpoint for cell-cycle entry); and when the signal is decreased to some critical value Σ 1 , the cell returns to its initial (post-mitotic) phase.We shall now describe our computational results for the full system (1)- (5).Any choice of initial values will become negligible after some time; for simplicity we take initial values Σ First, for the PRVs, we get the results in Figures 3 and 4. The curve spanned by dots in Figure 3 is the projection of the solution on the (C 2 , Σ) plane.Figure 4 shows the periodic behavior of all the components of the solution of (1)-( 5) after some time (approx.100 time units), so the effect of the choice of initial values becomes negligible.Note that A 2 < 1 2 , so the cell does not go into apoptosis.The nullcline of (1) is determined by the equation This Σ-nullcline varies with the parameters k s , k sd1 , and k sd2 .In particular, as we vary the signal input rate k s (and keep all the other parameters fixed), the Σ-nullcline changes, but the C 2 -nullcline does not.If, roughly, k s < 0.00065, and all other parameters are at PRVs then the two nullclines intersect below the lower knee of the S-shaped C 2 -nullcline, and the cell is quiescent (see Fig. 5).In this case the cell "gets stuck" at the intersection of the nullclines of Σ and C 2 before reaching the lower knee of the S-shaped curve.This is a quiescent cell, with A 2 < 1 2 (no apoptosis).If k s > 0.00251 then the two nullclines intersect above the higher knee of the Sshaped C 2 -nullcline (see Figure 6).The variable C 2 jumps up (i.e., the cell crosses the checkpoint for cell cycle entry), but it does not complete the cycle: it "gets stuck" at the intersection of the two nullclines.In this case A 2 > 1 2 after about 350 time units, i.e., the the cell enters apoptosis.For k s = 0.00065 the two nullclines intersect just above the lower knee of the C 2 -nullcline (see Fig. 7a), and for k s = 0.00251 they intersect just below the higher knee (see Fig. 8a).In both cases the system exhibits periodic oscillations.As shown in Figures 7b and 8b, for k s just above the quiescent limit, the system spends a lot of time in the inactive state (i.e., C 2 ≈ 0), and for k s just below the apoptotic limit the system spends a lot of time in the active state (i.e., C 2 ≈ 1).
As we can see in the bifurcation diagram in Figure 9, these oscillations are born out of Hopf bifurcations.Away from the point B, each vertical line intersects the near-rectangular thick curve at two points, the distance between them being the amplitude of a stable periodic solution.Near B, the thick curve actually bulges slightly leftward, so that each vertical line intersects it at four points.The distance between the two nearest points is the amplitude of an unstable periodic solution, and the distance between the two farther points is the amplitude of a stable periodic solution.
As shown in Figure 3, these oscillations are hysteretic in the sense that they involve periodic switching between the upper and lower branches of the C 2 -nullcline.Strictly speaking, these oscillations cannot be interpreted as the CDK oscillations in the cell cycle.However, it is interesting to note that the hysteretic oscillations, as illustrated in Figure 3 resemble the hysteretic oscillations of biological control systems which occur, for example, in the yeast cell cycle (see [9] for a review).For our purposes, we consider entry into the cell cycle when C 2 is switched on.
4. The role of other parameters.Consider the parameter k 2a , related to the input provided by the cell-cycle marker C 2 to the control node G 2 .This parameter measures the "depth" of the hysteretic loop of the C 2 -nullcline.If we increase k 2a from its reference value 0.01 to k 2a = 0.012, we get a "deeper" C 2 -nullcline, as we see by comparing the C 2 -nullclines from section 3 with the C 2 -nullcline in Figure 10a.On the other hand, if we decrease k 2a , then the "depth" of the S shape decreases.In particular, if k 2a = 0, then the C 2 -nullcline loses its S shape, so we cannot get a hysteretic loop (see Fig. 11).In this case the solution, instead of being periodic, converges to some fixed point at the intersection of the nullclines.
We shall now consider the case where k 2a = 0.012.If we keep all the other parameters at their reference values then, as seen in Figure 10a, the cycle begins its upward jump, but does not come down to completion, since the intersection between the two nullclines is now on the upper branch of the C 2 -nullcline.Therefore, the cell goes into apoptosis (A 2 > 1 2 ) after some time.However, if k s is decreased (from the reference value k s = 0.002 to k s = 0.001), then the cell goes into cycling (with A 2 < 1 2 ) as shown in Figure 10b.The above considerations suggest that the assertion A in section 2 is valid only for k 2a restricted to some interval k 2a < k 2a < k 2a , which includes the points 0.01 and 0.012.
If k 2a becomes very large then the two nullclines will always intersect along the stable upper branch of the S-shaped nullcline, and the trajectory will get stuck at that stable fixed point.In this case we can only have quiescence or apoptosis (i.e., there can be no cycling for any values of k s and k sd2 ).
Finally, let us look at the interplay between the parameters k sd1 , k sd2 , and k s .The three possible modes of the cell are summarized in the diagrams of  12a is determined by the condition that the Σ-nullcline passes exactly through the lower knee of the S-shaped curve (i.e., the lower knee of the C 2 -nullcline).Similarly, the equation of the higher line in Figure 12a is determined by the condition that the Σ-nullcline passes exactly through the higher knee of the C 2 -nullcline.If k sd1 is very small (e.g., k sd1 = 0.001, its reference value) then it is roughly the quotient k s /k sd2 that determines the mode (see Figure 12a).As we increase k sd1 both lines are raised, but the lower line is raised faster than the higher line.Hence, the cycling region shifts upward and at the same time decreases (see Figure 12b), and a new region appears between the two lines.For example, if all parameters are set at their reference values, except for k sd1 = 0.03, then the cell is quiescent as in Figure 5, instead of cycling as in Figures 8 and 9. 5. Conclusions.It has been shown experimentally (see [2,3,10]) that there is an optimum range of transcriptional activities for E2F-1 and c-Myc for cells to proliferate.If the intensity of growth factor signaling is too small, cells are quiescent and do not cycle.However, if signals are too large then the apoptosis program is triggered; this can be interpreted as an important fail-safe mechanism to prevent uncontrolled cell proliferation as it occurs in cancer.
A simple model developed in [1] established the above biological observation for a certain range of parameters.In this paper, we carried out a full analysis of the model for a wide range of parameters.
We have shown that the positive feedback loop (shown as "a" in Fig. 1, and represented by the parameter k 2a in equation ( 3)) is required for an S-shaped C 2nullcline.This nullcline geometry is essential for a switch-like behavior between the various cellular states, particularly from quiescence to cell cycling and/or apoptosis.The biochemical basis of the loop "a" is well established (see [1]); the transcriptional activities of the factors E2F-1 and c-Myc (represented by G 2 in the model) are increased by CDK2 (represented by C 2 in the model) via inhibition of the retinoblastoma protein (Rb) which inhibits the G 2 factors.Our analysis suggests that the parameters involved in this positive feedback loop are the key parameters for controlling the cellular state transitions.Our simulations on the effect of signal strengths are consistent with experimental reports on the induction of apoptosis (see [1]).
The cross-shaped phase diagram in Figure 12b is an interesting result.The diagram shows that the series of cellular state transitions depends on the magnitude of the negative feedback loop between C 2 and the signal Σ (as measured by the value of k sd2 ).For example, if the signal rate k s is increased while keeping k sd2 = 0.05, then the system goes from quiescence to cycling, and eventually to apoptosis.But for k sd2 below 0.02, say k sd2 = 0.01, the same increase in k s will not give rise to an intermediate cell-cycling state.The region denoted by "q-a" in Figure 12b is a region where the cell could either be quiescent or apoptotic depending on the cell's history.If the cell is previously quiescent and the signal increases to the "q-a" region, then the cell remains quiescent; however, a further increase in the signal will send the cell directly into apoptosis (without going through the intermediary cell cycle phase).Thus, we predict that: Assertion B. If the negative feedback loop between the cell cycle and signaling pathways that impinge on the E2F-1/c-Myc node is very weak (namely, k sd2 is small) then for large enough signal degradation rate (k sd1 ), the cell will bypass the cycling phase, as it goes (with increasing signal input rate k s ) from quiescent phase to apoptosis.
As we have mentioned in section 3, the periodic oscillations generated by the model may not correspond to experimentally observed CDK oscillations, but they resemble the oscillations of biological control systems which occur, for example, in the yeast cell cycle (see [9]).Using a different model, consisting also of five ODEs, Obeyesekere et al. [8] exhibited a bifurcation diagram similar to Figure 7, but for different variables; that model does not include apoptosis.
It would be interesting to extend our model to a more realistic network that includes key components of the cell-cycle module, such as Rb, cyclin E/CDK2, Cdc25A and p27Kip1 (see [1]); such an extension might enable one to test our assertion B experimentally.It would also be interesting to include in the model details of the apoptosis module that incorporate a positive feedback loop between executioner caspases and cytochrome c release, which was also reviewed in [1].Such a model would involve a larger system of ODEs.Cell 85:537-548.E-mail address: gcraciun@mbi.ohio-state.eduE-mail address: bdaguda@gmail.comE-mail address: afriedman@mbi.ohio-state.edu

Figure 1 .
Figure 1.Modular structure of the model.

Figure 2 .
Figure 2. Mechanistic implementation of the model.The solid arrows indicate input into or output from a node; dotted arrows indicate influence of a node on a process, without loss of mass.

Figure 3 .
Figure 3.The projection of the trajectory on the (Σ, C 2 ) plane at PRVs.

2 Figure 7 . 2 Figure 8 .
Figure 7.The case where k s is near the quiescent limit: (a) reduced phaseplane (b) explicit functions of time.

Figure 9 . 2 Figure 10
Figure 9. Bifurcation diagram of C 2 (vertical axis) with respect to k s (horizontal axis).The values of all other parameters are given by the PRVs.The Hopf bifurcation at A is supercritical, and the Hopf bifurcation at B is subcritical.The thin solid lines represent stable steady states, and the thin dashed line represents unstable steady states.

Figure 13 .
Figure 13.Pathways linking the transcription factors c-Myc and E2F-1 with the initiation of the cell cycle and apoptosis.Dashed lines are transcriptional processes; solid lines are post-translational processes.Arrows signify activation, while hammerheads signify inhibition.