Models of Spatially Restricted Biochemical Reaction Systems

Many reactions within the cell occur only in specific intracellular regions. Such local reaction networks give rise to microdo-mains of activated signaling components. The dynamics of microdomains can be visualized by live cell imaging. Computational models using partial differential equations provide mechanistic insights into the interacting factors that control microdomain dynamics. The mathematical models show that, for membrane-initiated signaling, the ratio of the surface area of theplasmamembranetothevolumeofthecytoplasm,thetopol-ogyofthesignalingnetwork,thenegativeregulators,andkinetic properties of key components together define microdomain dynamics. Thus, patterns of locally restricted signaling reaction systems can be considered an emergent property of the cell.

Many reactions within the cell occur only in specific intracellular regions. Such local reaction networks give rise to microdomains of activated signaling components. The dynamics of microdomains can be visualized by live cell imaging. Computational models using partial differential equations provide mechanistic insights into the interacting factors that control microdomain dynamics. The mathematical models show that, for membrane-initiated signaling, the ratio of the surface area of the plasma membrane to the volume of the cytoplasm, the topology of the signaling network, the negative regulators, and kinetic properties of key components together define microdomain dynamics. Thus, patterns of locally restricted signaling reaction systems can be considered an emergent property of the cell.
Homeostasis and regulated state change in the living cell require many processes to function in a coordinated fashion. Such multitasking is precisely orchestrated (1). For this, the cell integrates information from many extracellular signals in a timely and specific manner to evoke physiological processes. Complex networks of interacting pathways and regulatory mechanisms control and modulate information flow (2). Such complex systems cannot be intuitively understood. Mathematical modeling has served as a powerful analytical and predictive tool to understand how signals evoke responses and identify emergent properties of networks (3) that create an ensemble of output patterns that control cellular response machines. In this review, we describe how computational models help us understand the biochemical mechanisms underlying spatially restricted signaling events within cells.

Functional Specificity from Spatially Restricted Signaling
Organelles often spatially constrain cellular functions. Electrical activity occurs at the plasma membrane, ATP production in the mitochondria, and transcription in the nucleus. It is increasingly recognized that even when such organelle-based barriers do not exist, much of the biochemistry in the cell functions in a spatially restricted manner. The best examples of such spatially restricted biochemical reaction systems come from cell signaling pathways. Many receptors that evoke diverse cellular responses activate a relatively smaller number of shared intracellular signaling pathways, raising the question of how signaling specificity is achieved. One mechanism for achieving response specificity is the spatial separation of signaling reactions. The notion that signaling was spatially restricted was first proposed for the cAMP system by Brunton et al. (4) over 20 years ago. Compartmentalization was skeptically received then, given the popularity of the fluid mosaic membrane model (5) at the time and the prevalent idea that diffusible molecules would be impossible to segregate.
An example of spatially restricted signaling is found in cardiac myocytes, where elevation of cAMP levels regulates contractility, electrical excitability, and carbohydrate metabolism in a hormone-specific manner. For instance, ␤ 2 -adrenergic stimulation leads to an increase in myocyte contractibility due to the PKA 2 -dependent modulation of L-type Ca 2ϩ channels and troponin I (6), an actin-binding protein that regulates the myosin ATPase activity via tropomyosin (7). In contrast, other ligands such as prostaglandin E 1 can induce similar temporal increases in cAMP but are incapable of regulating cardiac muscle contraction. These data suggested that cells spatially and functionally compartmentalize intracellular signaling. Advances in live cell imaging, new fluorescent probes, and use of fluorescence resonance energy transfer for the visualization of biochemical reactions have shown that, in the live cell, spatially localized signaling is a common phenomenon. Using these approaches, it has been shown in several cell types that increases in cAMP levels are local (8,9).
Local signaling can provide context specificity to cellular responses to the MAPK pathway. Activation of Ras can lead to differentiation or proliferation of PC12 cells. How these two opposing cellular events can result from activation of the same pathway has been explained by correlating the response with duration of MAPK activation (10), which in turn depends on the location of Ras activation. Ras is activated at two sites within the cell: the plasma membrane and the Golgi membrane (11). Activation of Ras at the Golgi is prolonged and leads to differentiation, whereas brief plasma membrane activation leads to proliferation. The determinant of the location of Ras activation is the length of the intracellular Ca 2ϩ signal, as a prolonged Ca 2ϩ signal leads to translocation of the Ras guanine nucleotide exchange factor RasGRP1 to the Golgi membrane, where it activates Ras (11).
During maturation, the fate of thymocytes is determined either by their ability to recognize themselves (negative selection) or by the presence of a functional T-cell receptor (positive selection). Both processes are dependent on the Ras-MAPK pathway. Location of Ras activation determines whether a thy-* This work was supported, in whole or in part, by National Institutes of Health Grant GM 072853 and Systems Biology Center Grant P50 GM 065085. This is the fourth article of six in the Thematic Minireview Series on Computational Biochemistry: Systems Biology. This minireview will be reprinted in the 2009 Minireview Compendium, which will be available in January, 2010. □ S The on-line version of this article (available at http://www.jbc.org) contains supplemental Fig. S1. 1 To whom correspondence should be addressed. E-mail: ravi.iyengar@ mssm.edu. mocyte will undergo negative or positive selection (12). Positive selecting ligands induced the recruitment of RasGRP1 to the Golgi membrane, leading to Ras activation at the Golgi, whereas negative selecting ligands promoted Ras activation at the plasma membrane. Selective activation of Ras in the Golgi, but not the plasma membrane, has been linked to T-cell activation (13). Cdc42 activation is also spatially delimited (14). In a fibroblast spreading on a fibronectin-coated surface, Cdc42 regulates Arp2/3 nucleation of actin at the leading edge. Visualization of Cdc42 activation showed that the region of the cell perimeter undergoing rapid spreading had an enhanced accumulation of active Cdc42, whereas filopodia had little or no active Cdc42 (14). These observations led to the concept of "microdomains" defined as functional regions of regulated local increases in concentrations of activated signaling components.

Intracellular Microdomains
What is a microdomain, and how are its dynamics regulated? A microdomain is a region of micron/submicron dimensions within which increases in levels of signaling components are controlled by the relative rates of its production, diffusion, and destruction. The number of reactions underlying a microdomain can range from at least two, as in the simplest case of Ca 2ϩ , to tens, as in the case of MAPK microdomains. A microdomain has two defining characteristics: a point at which the concentration of the activated signaling component is the highest and a gradient that defines the change in concentration of the active component from the highest point to the level present in the surrounding milieu. The slope of the gradient defines the boundaries of the microdomains. Both these characteristics are dynamic, i.e. they are changing continually with time. Hence, all microdomain views are snapshots in time, and steady-state assumptions are mostly invalid.
The idea of signaling gradients was first articulated by Turing (15) more than 50 years ago in the context of developmental biology, where extracellular morphogen gradients provide positional cues for a developing tissue. These gradients can be produced by the local synthesis of a diffusible molecule and its slow degradation across space-creating gradients that may be several microns long. Mechanisms underlying morphogen gradients have been widely studied experimentally and computationally (16,17). Analysis of intracellular gradients lagged behind both theoretically and experimentally. Since its discovery (18,19), the use of green fluorescent protein (20) and similar probes to track proteins and interactions in real time within live cells has accelerated progress.
Among the first intracellular microdomains to be visualized were those for Ca 2ϩ at the intracellular terminus of Ca 2ϩ channels because the speed of conductance of Ca 2ϩ channels is several orders of magnitude higher than the diffusion coefficient of Ca 2ϩ (21). Ca 2ϩ /calmodulin feedback inhibitory mechanisms (22) tightly control the buildup of Ca 2ϩ within these microdomains. Dissipation of these microdomains is due not only to diffusion but also to multiple sequestration mechanisms involving the endoplasmic reticulum stores and mitochondria. Several other types of spatiotemporal patterns of intracellular calcium regulation have been described (23).
Diffusible second messengers (24,25), activated small GTPases (14, 26 -29), and protein kinases (8,9,30) exist within microdomains. Factors underlying the formation and dissipation of these microdomains are complex. An integration of experimental data and mathematical modeling is needed to provide a mechanistic understanding. Even though experimental approaches yield physical evidence for the existence of microdomains, they offer only limited insight into origins. Mathematical models can complement the empirical approaches by probing possible mechanistic bases of the observations. Simulation-derived hypotheses can then be validated experimentally. The combination of experiments and mathematical modeling has highlighted a number of factors that contribute to the formation of microdomains. These factor are cell shape, topology of the reaction network, diffusion coefficients, and reaction kinetics (intrinsic properties of the components involved). To further complicate matters, the factors are interdependent. Thus, for a mechanistic understanding, it is necessary to determine how interdependence among these factors controls microdomain dynamics.

Computational Approaches to Modeling Spatially Restricted Biochemical Reactions
Models of spatially restricted biochemical reactions can be built using ODE-based compartmental models or PDE-based spatial models. ODE models describe the rate of change in the amount of molecules over time, and PDE models describe the rate of change in the amount of a molecule over time and space (36). Both approaches yield predictive models. Compartmental models can be easily solved using commercial simulators such as MATLAB. Modeling and running numerical simulations of realistic spatial models are more challenging, and there are very few simulators that biological researchers can use. Virtual Cell (37) is a freely available simulator useful for developing simulations from cell imaging data using both ODE and PDE models.
Compartmental models disregard details of spatial organization and the movement of molecules. The steps in setting up a compartmental model are shown in supplemental Fig. S1. Spatial information is abstracted to two entities, volume and surface area, and within each compartment, all molecules are assumed to be evenly distributed (i.e. fast diffusion). Even with these simplifying assumptions, compartmental models provide valuable insight into complex phenomena. A multicompartment model of myocyte function integrating mitochondrial biochemical reactions with electrical activity (ion channel function) in the plasma membrane and excitation-contraction in the sarcolemma provided valuable insight into how feedback mechanisms between mitochondrial oxidative phosphorylation reactions (cellular biochemistry) and the electrical and mechanical activities (cell physiology) drive metabolic regulation of cardiac function (38).
Spatial models deal with cell shape explicitly. To understand the regulatory insights spatial models provide, it is useful to know how these models are developed and computed. The steps involved in setting up a model with a realistic geometry in Virtual Cell are shown in supplemental Fig. S1. The reactions are set up as an ODE compartmental model, and then compartments are mapped onto a spatial geometry from a microscopic image. Mapping involves placing the components in the appropriate location at the appropriate concentration and defining their diffusion coefficients. Numerical simulations are then carried out using a finite volume regular grid solver, and results can be displayed as color-coded concentrations within the geometry, allowing for a visually accessible comparison between experiments and simulations.

Insights from Spatial Models
Spatial PDE models of cellular reactions are few because the diffusion coefficients of few cellular components are known. Although diffusion coefficients can be calculated from molecular weights, it is likely that molecular crowding will affect diffusion within the cell (39 -41). With fluorescence recovery after photobleaching experiments, it is becoming feasible to obtain direct measures of diffusion of cellular components. Systematic variation of the diffusion parameter can be used to identify computationally parameters that are likely to yield experimentally observed patterns of microdomain dynamics. Despite these limitations, spatial models have yielded the first mechanistic insights. Spatial models have shown that microdomain dynamics is a systems level property arising from the interplay of several variables. Mathematical models have identified spatial information as an entity that is transmitted through a signaling network separately from information regarding activity states. Here, we describe spatial models of microdomains in myocytes and neurons to highlight the insight obtained from these models.
An experimentally observed cAMP microdomain from the studies of Zaccolo and Pozzan (9) is shown in Fig. 1A. Results from a spatial model that simulate these observations are shown Fig. 1B. To understand what a computational model tells us, consider the schematic in Fig. 1C. A microdomain depends on the localization of a key upstream component, in this case, adenylyl cyclase, which is an intrinsic membrane protein.
Whereas cAMP-degrading phosphodiesterases are located largely in the cytoplasm, some phosphodiesterase isoforms are anchored to the endoplasmic reticulum or other membranes through anchoring proteins. The spatial separation of cAMP production and degradation results in a microdomain. The spatial gradient of cAMP is shown in Fig. 1D. In the myocyte, there is an extensive, specialized endoplasmic reticulum that forms a physical barrier for cAMP diffusion and thus defines the observed shape of the microdomains. Although contributing to microdomain characteristics, such a barrier is not essential.
Many protein kinases and phosphatases have anchors and scaffolds that create local signaling hubs. Prominent among these are A-kinase anchoring protein for PKA, microtubule affinity-regulating kinases for protein kinases C, JNK-interacting proteins for JNK, and MP1 for MAPK. Microdomains within the plasma membrane such as lipid rafts and caveolae (31) also promote spatial organization and separate extracellular signals spatially before their transduction. Certain signaling proteins such as lipid-linked non-receptor tyrosine kinases appear to be preferentially localized to these intramembrane structures (32). Intuitively, one would assume that such anchors are essential for microdomains; however, mathematical analyses indicate they are not essential.
When the signal originates at the membrane and is degraded in the cytoplasm, cell shape or the ratio of the surface area of the plasma membrane to the cytoplasmic volume (surface/volume ratio (S/V)) becomes important. The S/V of thin protrusions like filapodia or dendrites is larger than that of thicker structures such as the cell body of a neuron. Models predict and experiments show that, in regions of high S/V ratio, elevation of cAMP is enhanced (33), as is activation of Cdc42 and protein kinase C (34,35). Because both amplitude and duration of a signal are related to the numbers (moles) of signaling molecules along a pathway, increasing the surface area of the plasma membrane relative to the volume of the cytoplasm decreases the volume within which the signaling molecule can diffuse and the gradient dissipate. In most cases, as the concentration (molecules/volume) stays constant, increasing the S/V ratio increases the signal source. This results in a buildup of the activated signaling component and the resultant microdomain. Although computational analyses indicate that anchors and scaffolds are not required for the formation of microdomains, the role of signaling complexes in microdomain formation and dynamics has not yet been extensively studied. An intriguing hypothesis is that cytoplasmic anchors may help override some  Zaccolo and Pozzan (9) reprinted with permission from AAAS. The fluorescence resonance energy transfer ratio images are color-coded, with green signifying basal concentrations of cAMP and yellow to red signifying high concentrations of cAMP. Scale bar ϭ 10 m. The white arrow in the enlarged image points to a cAMP microdomain (ϳ1 m). B, simulation in Virtual Cell of receptor-triggered cAMP increases using the cell shape traced from the myocyte image in A. The black arrow shows the corresponding cAMP microdomain in same region as in A. C, schematic diagram of factors controlling cAMP microdomains in myocytes: plasma membrane (PM) localization of adenylyl cyclase (AC), separation of production and degradation (phosphodiesterase (PDE)) activities for cAMP, and diffusion hindrance by intracellular organelles such as the sarcoplasmic reticulum (SR, depicted as the black double dashed line). D, microdomain characteristics are defined by the highest concentration point within the microdomain and the resulting gradient that extends to the periphery and defines the microdomain boundary. These are represented as a function of distance from the site of production. Arbitrary numbers are used for illustrative purposes.
of the spatial constraints that low S/V areas impose on microdomains. This hypothesis can be computationally addressed.
Computational analyses show that S/V ratio is only one of several factors regulating microdomains. Negative regulators also play a critical role. Dependence of cAMP microdomains on the activity of phosphodiesterases has been shown (9) as depicted in Fig. 1C. At its core, this phenomenon is analogous to Turing's original description of a point source with a diffusing signal and negative regulator resulting in a gradient, and a signaling version of this phenomenon was described (42). The toy graph in Fig. 1D shows this gradient. From S/V ratio consideration, the notion arises that thin, elongated regions of a cell (such as filopodia or dendrites in neurons) are favored for the formation of microdomains. This is shown in the simulations in Fig. 2. Here, a cAMP-PKA-MAPK network in a neuronal den-dritic arbor was modeled ( Fig. 2A). Although the receptor, G s , and adenylyl cyclase are uniformly distributed throughout the cell, the thin dendrites (high S/V ratio) show the highest levels of PKA activation (Fig. 2B). Thus, cell shape becomes a determinant of microdomain characteristics (33,34). When the phosphodiesterase is fully inhibited, then irrespective of cell shape, PKA activation is uniformly elevated (Fig. 2C). These observations indicate that the biochemistry of the system can override the physical constraints.
Although S/V ratios, cell shape, and negative regulators explain the dynamics of microdomains of signaling molecules such as cAMP, Ca 2ϩ , and even some small GTPases (34), these two factors do not explain the gradients in downstream components such as MAPKs. The computational analyses indicate that network topology is also important (Fig. 2D). The hierarchy of the negative regulators within the network topology determines microdomain characteristics of the downstream component. Both the upstream cAMP phosphodiesterase and downstream PTP control the MAPK microdomain. This control is evident only when S/V conditions are high; thus, cell shape becomes an overall constraint within which the network topology can be altered to determine downstream microdomains. These computational studies showed that the microdomain characteristics (i.e. gradient) of PKA are transmitted to MAPK through PTP but not through b-Raf and MEK. Nevertheless, the magnitude of the increase in MAPK activity results from b-Raf and MEK activation (Fig. 2D). Thus, information regarding position of the activated component (i.e. spatial information) is handled differently from information regarding activity state. These two types of information can be transmitted through different paths within a signaling network (Fig. 2D). This insight, originally obtained from the computational model, was experimentally verified by ablating the key component (PTP) utilized for spatial information transfer. When PTP was ablated, information regarding activation state was still transmitted, but the spatially restricted activation of MAPK was lost, indicating that information regarding activity state is distinct from spatial information (33).
Separation of spatial information from activity information raises an interesting question: what are the characteristics , PKA microdomains when phosphodiesterase is active. D, network topology and kinetic parameters control the transfer of microdomain characteristics (i.e. spatial information) from upstream to downstream components. These simulations show that kinetic parameters are crucial factors in determining the transfer of microdomain characteristics. In the model, PKA has a high K cat for b-Raf, which leads to extensive activation of b-Raf and loss of microdomains. In contrast, PKA has a low K cat for PTPSL, resulting in microdomains of phosphorylated PTPSL (PTPSL i ) that are similar to those of PKA. The regulation of PTP by PKA* leads to two populations of phosphatase activity, PKA*-inhibited activity (PTPSL i ) and the active non-phosphorylated enzyme (PTPSL*). Areas with a high concentration of PKA* have low PTP activity and consequently high phospho-MAPK. Areas with a low concentration of PKA* have high PTP activity and low phospho-MAPK. This spatial and activity heterogeneity of PTPSL results in recapitulating the PKA* microdomains as MAPK* microdomains. within a pathway that allow for the flow of spatial information? Mathematical simulations reveal that the relationship between the kinetic parameters of the upstream and downstream components is an important determinant. If the upstream enzyme extensively activates (high K cat ) the downstream component (e.g. PKA 3 b-Raf) (Fig. 2D), the microdomain characteristics of the upstream component will not be preserved in the downstream component. Inversely, a low K cat reaction such as the phosphorylation of PTP by PKA results in the phosphorylated form of PTP exhibiting microdomains similar to those of upstream PKA* (Fig. 2D). The non-homogenous distribution of PKA* results in differing ratios of active (non-phosphorylated PTPSL*) to inactive (phosphorylated PTPSL i ), allowing for the formation of activated MAPK1/2 microdomains (Fig. 2D). Thus, the microdomains of PKA* are recapitulated in the microdomains of phospho-MAPK. The extent to which the signal controls the negative regulator becomes a critical factor. Too much active PTP in all locations would result in a nearcomplete loss of phospho-MAPK, and too much PTP inhibition would result in high levels of phospho-MAPK and in the dissipation of the microdomain. This is a Goldilocks relationship in which the kinetic parameters and network topology need to be tuned just right for the spatial information to be transmitted. The balance between kinetic parameters and network topology works only if the overall S/V ratios are in the correct range. For spatial information to flow through the cell, the cell shape, signaling network topology, and kinetic characteristics of the signaling components need to be balanced. Hence, spatial information flow is an emergent property that arises from the functional organization of the cell within an optimized shape.

Perspective
Many questions regarding the role of scaffolds and anchors in the dynamics of microdomains remain to be addressed. The study by Fuller et al. (43) on the spatiotemporal dynamics of the Aurora B kinase phosphorylation of its substrates during anaphase highlights some of these questions. Interactions with anchors and scaffolds are often regulated by signaling. This regulation introduces the possibility that dynamic scaffolding can play a substantial role in controlling the flow of spatial information. This is an area for future studies where computation and experimentation can be effectively integrated.