Computational approaches for modeling regulatory cellular networks

Abstract

Cellular components interact with each other to form networks that process information and evoke biological responses. A deep understanding of the behavior of these networks requires the development and analysis of mathematical models. In this article, different types of mathematical representations for modeling signaling networks are described, and the advantages and disadvantages of each type are discussed. Two experimentally well-studied signaling networks are then used as examples to illustrate the insight that could be gained through modeling. Finally, the modeling approach is expanded to describe how signaling networks might regulate cellular machines and evoke phenotypic behaviors.

Figure 1

Switching behavior (bistability) resulting from interacting signaling pathways. (a) By binding to α1 and α2 adrenoceptors, norepinephrine can activate both the G_q and G_i/o pathways. Interactions between these pathways can result in a positive-feedback loop [78]. A key question regarding the behavior of this system is whether the positive-feedback loop can lead to bistability. (b) To answer this question, we can develop a model of this system to determine whether bistability occurs. The first two ordinary differential equations (ODEs) describe the formation of receptor–ligand complexes for α1 and α2 adrenoceptors. The equations state that the rates of change of the receptor–ligand complexes are equal to the rate of binding of the ligand to the receptor (first term), minus the rate of dissociation of the ligand (second term), minus the rate of binding of the G proteins to the complexes, plus the rate of activation of the G proteins. In these equations, the rate of change in concentration is represented as d[C]/dt, where [C] represents the concentration of species C. The parameter k_i represents the kinetic constant (forward or reverse rate) for reaction i. The last equation of this set states that the rate of change of phosphorylated mitogen-activated protein kinase (MAPK) is equal to the rate of phosphorylation by MEK minus the rate of dephosphorylation. Here, the phosphorylation of MAPK is assumed to follow Michaelis–Menten kinetics, with the rate specified by the affinity constant K_m and the turnover number k_cat of MEK. Similar equations can be written for the remaining components of the system. After the equations have been written for all the components, and values assigned for concentration and rate terms, the system can be analyzed for many properties. For instance, bistability can be explored by examining the effect of kinetic parameters and initial concentrations on the steady-state solutions, which are calculated by setting the rate of change of all molecules to zero. Abbreviations: AA, arachidonic acid; AR, adrenoceptor; cPLA₂, cytoplasmic phospholipase A₂; DAG, diacylglycerol; IP₃, inositol (1,4,5)-trisphosphate; NE, norepinephrine; PKC, protein kinase C; PLCβ, phospholipase Cβ. Image in (a) reproduced, with permission, from Ref. [78].

Figure 2

Modeling spatially restricted biochemical reactions. (a) Epidermal growth factor (EGF) activates the canonical Ras mitogen-activated protein kinase (MAPK) pathway on the plasma membrane through Grb2/Sos. In addition, EGF can trigger an elevation in intracellular Ca²⁺ through phospholipase Cγ (PLCγ). The exchange factor RasGRP1 translocates to the Golgi after Ca²⁺ binding, activating Golgi-bound Ras. Intracellular Ca²⁺ also activates the GTPase-activating protein (GAP) CAPRI, which translocates to the plasma membrane and deactivates plasma-membrane-bound Ras. A key question in this system is whether Ras molecules on the plasma membrane and those on the Golgi are activated for different durations by the EGF receptor (EGFR). (b) To answer this question, we need to build a model that includes spatial specifications. Hence, we use a set of partial differential equations (PDEs). The first two equations in this figure describe the activation in the cytoplasm, diffusion and capture of RasGRP1 on the Golgi membrane. The first equation describes the rate of change of activated RasGRP1 with respect to time. The rate of change of activated RasGRP1 has three components. The first term describes the diffusion of RasGRP in the cytoplasm, which is proportional to the second derivative of the concentration, (1/r)∂/∂r(r∂[Ca_RasGRP*]/∂r), multiplied by the diffusion coefficient D_RasGRP*. The second and third terms describe the activation of RasGRP by binding of Ca²⁺ and deactivation by release of bound Ca²⁺. The parameters k₁ and k₂ are the rate constants for these processes. The second equation is called a boundary condition; in this case, it states that the flux of activated RasGRP1 to the Golgi membrane, 2πR∂[Ca_RasGRP*]/∂r, is equal to the rate of binding of RasGRP to Ras, k_onCa_[RasGRP*]|r=R, where the Golgi is defined as a circle with a radius (r) of R. The last equation states that the rate of change of RasGTP on the Golgi is equal to the rate of activation of Ras by RasGRP minus the rate of inactivation by a GAP. Note that, for clarity, the diffusion of RasGTP in the plane of the Golgi membrane is ignored. Similar equations and boundary conditions must be specified for the rest of the components. Analytical solutions are not guaranteed for many systems of PDEs, and PDE models are typically solved numerically. Abbreviation: IP₃, inositol (1,4,5)-trisphosphate.

Figure 3

Feed-forward relay for decoding signal duration. Recent studies by Blenis and co-workers have elucidated a potential mechanism for decoding the duration of a mitogen-activated protein kinase (MAPK) response [53,54]. The decoder essentially works as a feed-forward relay system. Activated MAPK induces transcription of intermediate-early genes such as c-fos. After c-fos is transcribed and translated, it can be degraded by the ubiquitin–proteasome pathway or stabilized by MAPK-1,2 interaction and phosphorylation. If the duration of MAPK activation is short compared with the lag associated with c-fos transcription and translation, there will be little activated (MAPK*) to phosphorylate the newly synthesized c-Fos protein, and c-Fos is ubiquitinated (c-Fos–Ub) and degraded. However, if the duration of MAPK activation is longer than the time required for transcription and translation of c-fos, MAPK* can bind to and phosphorylate c-Fos (c-Fos–P) and protect it against degradation. In turn, the stable c-Fos stimulates the transcription of Fos-dependent genes, even after MAPK is deactivated and, thus, the initial signal is propagated across time scales.

Figure 4

Modular organization of regulated cell movement. A major goal of mathematical modeling in cell biology is to predict cellular function from biochemical reaction networks. Such models might involve multiple networks and can be studied best in a modular fashion. Regulated cell movement can be described by two interacting modules comprising the regulatory network (green background) and the motility machines module (yellow background), which consists of the cytoskeleton and molecular motors that, together, make up the force-generating machinery. The signal network processes and integrates information from extracellular signals, in addition to internal signals from the force-generating machinery (mechanotransduction) to regulate the output force that leads to cell movement. The reverse arrows indicate the ability of cellular localization to regulate the capability of the signaling network, and the ability of force to affect directly both the function of the cytoskeleton and the motors (mechanochemical regulation) and to affect indirectly the signaling network (mechanotransduction).