Modeling the dynamic behavior of biochemical regulatory networks

Abstract

Strategies for modeling the complex dynamical behavior of gene/protein regulatory networks have evolved over the last 50 years as both the knowledge of these molecular control systems and the power of computing resources have increased. Here, we review a number of common modeling approaches, including Boolean (logical) models, systems of piecewise-linear or fully non-linear ordinary differential equations, and stochastic models (including hybrid deterministic/stochastic approaches). We discuss the pro's and con's of each approach, to help novice modelers choose a modeling strategy suitable to their problem, based on the type and bounty of available experimental information. We illustrate different modeling strategies in terms of some abstract network motifs, and in the specific context of cell cycle regulation.

Keywords: Bifurcation theory; Dynamic models; Logical models; Molecular regulatory networks; Piecewise-linear odes; Signaling motifs; Stochastic models.

Figure 1.

Regulatory network diagrams. (A) Five representative topological motifs. X, Y and Z are proteins. A barbed connector represents the ‘activation’ of one protein by another (perhaps by itself). A blunt connector represents ‘inhibition’. (B) Four possible biochemical reaction mechanisms that implement motif (i) in panel A. In the first case, X and Y inhibit each other’s production. In the second case, Y is a protease that degrades X, whereas X is a kinase that phosphorylates and inactivates Y. In the third case, two kinases phosphorylate and inactive each other. In the fourth case, X inhibits the production of Y and Y activates the degradation of X.

Figure 2.

Phase-plane portraits of piece-wise linear (PWL) and smooth, nonlinear (soft-H) ODE models of motifs discussed in the text. (A) PWL model of Motif 1A(i)-Case A, with γ₁ = γ₂. (B) PWL model of Motif 1A(i)-Case C, with γ₂ = 2γ₁. (C) Soft-H model of Motif 1A(i)-Case A, with γ₁ = γ₂. (D) Soft-H model of Motif 1A(i)-Case C, with γ₂ = 2γ₁. Parameter values are given in Table 3. In all four panels, the black lines are ‘trajectories’ starting from a variety of initial locations in state space and proceeding to one or the other of two stable steady states, at C₁=1, C₂=0 or at C₁=0, C₂=1. In cases C and D, the red curve is the C₁-nullcline (where dC₁/dt = 0) and the green curve is the C₂-nullcline (where dC₂/dt = 0). The red and green curves intersect at the third steady state, an unstable saddle point, in the interior of the unit square. The curve passing through the saddle point (yellow in panel C) separates the ‘domains of attraction’ of the two stable steady states.

Figure 3.

Phase-plane portraits (A, B) and one-parameter bifurcation diagrams (C, D) of soft-H ODE models of Motifs 1A(ii) (panels A and C) and 1A(iii) (panels B and D). Parameter values are given in Table 3; (A) S= 0.5, (B) S= 1, (C and D) S is the bifurcation parameter. In panel C we plot the steady state value of C₁ as a function of S. At S = 0, there is a single steady state (the ‘naïve’ state) at C₁ = C₂ ≈ 0. At S = 0.04, the system undergoes a pair of saddle-node bifurcations that create two new, ‘differentiated’ states (one at C₁ ≈ 0.82, C₂ ≈ 0 and the other at C₁ ≈ 0, C₂ ≈ 0.82) in addition to the naïve state. At S = 0.28, the naïve state is lost at a pitchfork bifurcation. The two differentiated states persist for 0.04 < S < 0.67; however, at S = 0.43, a second pitchfork bifurcation creates a new, stable steady state (called a ‘dual expressing’ state) for which C₁ = C₂ > 0.88. In panel D the system executes stable limit cycle oscillations for signal strengths between the Hopf bifurcation at S = 0.825 and a saddle-node-on-an-invariant-circle (SNIC) bifurcation at S = 1.30. There is a second saddle-node bifurcation at S = 3.64, but the node is unstable so the bifurcation is unobservable in experiments. (For a brief explanation of the terminology of bifurcation points, see Table 4.)

Figure 4.

Standard component model of mitotic control in a ‘generic’ eukaryotic cell. (A) Generic model. CycB = cyclin B-dependent kinase, which drives cell into mitosis; Cdh1 = active form of the cyclosome, which degrades cyclin B and drives cell out of mitosis; CAP = counter-acting phosphatase, which activates Cdh1 during exit from mitosis and entry into G1 phase; INH = a stoichiometric inhibitor of CAP. (B) XPP-ode file for an SCM implementation of the generic model in panel A. CycB is governed by an ODE like Eq. (12); Cdh1 and INH are governed by ODEs like Eq. (10); and CAP is given by a tight-binding approximation like Eq. (17). Progression through the cell cycle (G1-S-G2-M) is driven by growth. The variable M(t) is cell mass, which increases exponentially at a specific growth rate ‘mu’. The mass-doubling time for a cell is 0.693/mu. The cell divides (M → M/2) when CycB level drops below a threshold, theta. (C) Simulation of the SCM for the parameter values given in panel B. The simulated cell undergoes periodic divisions every 58 time units, which is exactly the mass-doubling time. The simulated cell persists in G1 phase for ~38 time units, and then proceeds through S-G2-M in ~20 time units. For a mammalian cell, we might take 1 time unit ≈20 min; for a yeast cell, we might take 1 time unit ≈2 min. (D) One-parameter bifurcation diagram (schematic) for the SCM, treating ‘mass’ as a parameter instead of a variable. Saddle-node bifurcations at M = 0.6055 and 1.3355 bound a region of bistability. Limit cycle oscillations are evident between a saddle-loop bifurcation at M ≈ 1.3 and a Hopf bifurcation at M ≈ 7.4. (E) Stochastic simulation of the model studied by a deterministic SCM in panel C. The turquoise curves show stochastic fluctuations in Cdh1_total and INH_total; the orange curve, CAPT. The red curve is CycB abundance. The black curve shows the fluctuating activity of Cdh1 (high in G1 phase, low in S-G2-M), and the blue curve shows how CAP is released from binding to INH as the cell exits mitosis. The saw-tooth black curve is 1000×mass. Notice that cell size at division and cell cycle time are stochastic properties of cell cycle progression in this model.