Bioinformatics Advance Access originally published online on October 11, 2006
Bioinformatics 2006 22(23):2918-2925; doi:10.1093/bioinformatics/btl497
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Cell++—simulating biochemical pathways
1 Hospital for Sick Children, 555 University Avenue Toronto, ON, M5G 1X8, Canada
2 Department of Molecular and Medical Genetics, University of Toronto Toronto, ON, Canada
3 Department of Biochemistry, University of Toronto ON, Canada
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT 1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSIONREFERENCES |
|---|
Motivation: With the generation of a wealth of information, detailing cellular components, their functions and interactions, there is a growing need for the development of new computational tools capable of interpreting these data within spatial and dynamic contexts. Here, we introduce Cell++, a novel stochastic simulation environment with the capacity to study a wide variety of biochemical processes within a spatial context.
Results: Focusing on three case studies, we highlight the potential impact of spatial organization in the evolution and engineering of signaling and metabolic pathways. In addition to altering signaling and metabolic efficiency, simulations also demonstrated features consistent with the phenomenon of metabolic channeling.
Availability: Cell++ is licensed under the GNU general public license (GPL) and has been successfully implemented under Linux and IRIX operating systems. Source code together with a simple tutorial is available at http://www.compsysbio.org/CellSim/.
Contact: jparkin{at}sickkids.ca
Supplementary information: Supplementary data for this paper are available at Bioinformatics online.
|
1 INTRODUCTION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSIONREFERENCES |
|---|
The advent of genomics and post-genomics technologies is leading to the generation of vast amounts of biological data from catalogs of genes, their products and expression profiles, to their interactions with other molecules, kinetics and cellular locations. Recently, attempts have been made to place these data within the context of higher order networks (Hermjakob et al., 2004; Salwinski et al., 2004). However, it must be appreciated that these networks are complex dynamical systems and understanding the relationships between the components is often non-intuitive. Therefore, to study systems, such as signaling and metabolic pathways, involving 10's, 100's or even 1000's of different species of components, requires the development and application of novel computational methods.
Previous attempts to model biochemical pathways have tended to rely on the use of ordinary and partial differential equations (ODEs and PDEs, respectively) (Tomita et al., 1999; Kierzek, 2002). However, such systems of equations are often difficult to solve and often neglect potentially important spatial and organizational factors that can significantly affect pathway behavior (Ingber et al., 1994; Ingber, 2003; Vermassen et al., 2003; Ovadi and Saks, 2004; Takahashi et al., 2005). For example, it has previously been suggested that effective signal transduction in mitogen activated protein kinase (MAPK) cascades cannot be explained by diffusion alone (Kholodenko, 2003), but may instead be facilitated through active molecular transport and component co-localization. Similarly, the co-localization of common enzymes from metabolic pathways leads to the phenomenon of ‘metabolic channeling’, increasing catalytic efficiency and preventing the build up of substrate intermediaries (Jorgensen et al., 2005). In addition to complex formation, such co-localization of pathway components may arise from the creation of compartments created through the orchestrated interplay of membrane and/or cytoskeletal elements (Verkman, 2002; Ovadi and Saks, 2004; Takahashi et al., 2005; Monastyrskaya et al., 2005).
The recent recognition of the importance of including spatial factors in simulations of biochemical systems, has led to the development of several new modeling applications including SmartCell (Ander et al., 2004), project CyberCell (Broderick et al., 2004), MesoRD (Hattne et al., 2005) and MCell (Coggan et al., 2005). Each has successfully been applied to the study of particular classes of systems, from simple signaling systems to more complex systems, such as the release of neurotransmitters in synapses and membrane growth. However, while these simulators may be adapted to explore alternative pathway architectures and spatial organizations, their ability to consider other types of biochemical systems, such as metabolism is relatively limited (Takahasi et al., 2005).
Here we present Cell++, a unique simulation environment that provides the user with the flexibility to model a variety of biochemical pathways from signal transduction to metabolism. Combining a cellular automata engine with Brownian dynamics, Cell++ is able to simulate the bulk properties of large quantities of small molecules (e.g. pyruvate and Ca2+), while simultaneously allowing more accurate representations of large molecules that can display more complex behavior. To demonstrate its capabilities, we present three case studies in which Cell++ is applied to examine the role of spatial organization in signaling and metabolic pathways. Results from these investigations reveal the potential impact of spatial organization, both in terms of pathway efficiency and in governing pathway behavior.
|
2 METHODS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSIONREFERENCES |
|---|
2.1 Implementation
Minimum system requirements for Cell++ are 1 GHz Pentium or equivalent with 256 MB RAM and 1 GB HDD running Red Hat Linux 9.0. It has also been tested on Fedora Core 3.0 and 4.0; Ubuntu (both Linux); and SGI IRIX 6.5. Simulations are visualized by the use of the OpenGL graphics library (http://www.opengl.org). The basic design of Cell++ features a hybrid Brownian dynamics/cellular automata approach in which discrete components (typically proteins) move within a continuum superimposed on an underlying 3D cubic lattice which describes localized spatial components including the cellular environment (e.g. membrane and cytosol) and the concentration of small molecules (Fig. 1). The size of the lattice is user defined and determines the spatial resolution of small metabolite concentrations. A fine resolution, e.g. 5 nm allows for a more accurate spatial representation of the simulation environment, including small molecule concentrations, but will require a greater number of lattice units for a given volume. On the other hand, a coarse resolution, e.g. 200 nm, can be used to reduce the number of lattice units increasing simulation speed, important for systems involving small molecule diffusion. Simulation time scales linearly with the number of lattice cells considered.
|
Within the simulation environment, different cell shapes and intracellular geometries are readily definable. Movement of discrete components is handled as random walks that may be restricted by the local environment (e.g. membrane boundaries). Discrete component interactions (e.g. phosphorylation events) are modeled using a Monte Carlo approach in which the probability of interaction depends upon the type of components, their relative distance and local environment. Small molecule diffusion is handled by deterministic sets of equations in which a proportion of each molecule (determined by its diffusion coefficient) in each lattice site is moved to adjacent lattice sites. Additional sets of equations describe changes in concentration according to interactions with other cellular components (e.g. enzyme-substrate turnover or active transport). Simulations proceed via a series of iterations (representing discrete time steps) in which component movement and interactions are updated. System specific behaviors, such as small molecule uptake or protein generation/turnover are readily incorporated by e.g. allowing a probability of such actions during each iteration. For further details of Cell++ implementation see Supplementary material.
Cell++ is a standalone package written in C++. A simple graphical user interface allows the view of the system to be manipulated and provides control over its progress. New simulation systems are created by modification of the source code and its corresponding input files. Two tutorials detailing how user defined signaling and metabolic pathways may be simulated are available with the source code at http://www.compsysbio.org. Although outside the scope of this current work, we have also implemented more complicated behavior, such as the formation of complexes and regulation of enzyme activities, the details of which we would be happy to share with interested users.
2.2 Details of studied systems
Full details of implementation for the three presented case studies: a simple signal transduction pathway, part of glycolysis and calcium-induced calcium-release (CICR) are provided in the Supplementary material.
2.2.1 Signal transduction.
We considered a system with five distinct types of signaling molecules: membrane receptors, three cytosolic components representing consecutive elements in a kinase cascade (kinases a-c) and nuclear pores. Each molecule was treated as a discrete component. Signal transduction was modeled as a series of ‘activation’ events (representing e.g. protein phosphorylation). Starting with the membrane receptors (initially all set as activated), each component is able to activate the following component in the pathway until the final component (kinase c) encounters a nuclear pore, after which the signal transduction is considered complete and the simulation terminated. A 30 x 30 x 30 cubic lattice was used to describe four environments: an outer ‘membrane’ (one lattice cell thick); an inner ‘nucleus’ composed of 8 x 8 x 8 lattice units; a ‘nuclear membrane’ (one lattice cell thick) surrounding the nucleus and a ‘cytosolic’ compartment located between the two membranes. For some simulations, cytosolic components (kinases a-c) were constrained within a box of dimensions 8 x 10 x 10 lattice units joining the nuclear and outer membranes.
2.2.2 Glycolysis (part of).
This system was composed of both discrete components—representing four types of enzymes: phosphoglycerate kinase; phosphoglycerate mutase; enolase and pyruvate kinase—and concentrations assigned to each point on a 10 x 10 x 10 cubic lattice—representing five substrates: 1,3-bisphosphoglycerate, 3-phosphoglycerate, 2-phosphoglycerate, phosphoenolpyruvate and pyruvate enzymes. The size of each lattice unit was defined as 50 nm. In certain simulations, we defined a ‘constraining’ environment at the center of this environment consisting of a single lattice site in which certain types of enzyme were restricted. To reduce system complexity, other small molecules involved in reaction processes (e.g. ATP) were considered to be in abundance. Enzyme kinetics were calculated with reference to one substrate-one product reversible Michaelis–Menten equations. The parameters used in these simulations are listed in Supplementary Table 1 (Teusink et al., 2000).
2.2.3 CICR.
Adopting a 40 x 40 x 10 cubic lattice, one half (40 x 40 x 5 lattice units) was divided into 16 equal compartments (intraluminal calcium stores), each of which consisted of 9 x 9 x 5 lattice units surrounded by a Ca2+ impermeable membrane 1 lattice cell thick. The remaining space (40 x 40 x 4 lattice units) was defined as cytosol, through which Ca2+ can freely diffuse. Ca2+-release channels were placed on the membrane separating the intraluminal stores from the cytosol in one of three spatial distributions: random, clustered and uniform (Supplementary material). Ca2+-release channels were modeled as discrete entities, with 10–250 used in each simulation. Ca2+ was modeled as a concentration at each lattice site. Ca2+-release channels did not move. Once triggered by a sufficiently high cytosolic Ca2+ concentration (2.5 U/lattice cell), the release channels allow intraluminal Ca2+ to diffuse across the membrane for a limited number of iterations (150). For each simulation, a CICR was initiated by placing an ‘open’ Ca2+ channel in one corner of the membrane. Simulations terminated once the concentration in the lattice cell in the opposite corner increased to 2.5 U/lattice cell or after 8000 iterations (after which CICR was assumed to have failed).
|
3 RESULTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSIONREFERENCES |
|---|
3.1 Case study 1—a signal transduction pathway
Signal transduction pathways, such as MAPK cascades, operate as networks of interacting proteins in which a signal may be transferred through the pathway by modification e.g. phosphorylation and subsequent spatial relocation (Kondoh et al., 2005). Previously it has been suggested that effective signal transduction, cannot rely on simple diffusion alone but may require the formation of complexes and/or active transport processes (Kholodenko, 2003). Further speculation suggests that pathway behavior may additionally be modified by the organization of the high density of proteins within the cell, commonly referred to as molecular crowding (Elowitz et al., 1999; Ovadi and Saks, 2004; Takahashi et al., 2005).
Here we model a simple MAPK based cascade to:
- explore the hypothesis that co-localization, leading to efficient signal transduction, may be achieved through sub-cellular compartmentalization and
- investigate the influence of alternative molecular crowding architectures.
As an initial trial of Cell++, we began exploring the effects of altering relative component concentrations. Figure 2A shows how the time of signal transduction, (which we define as the amount of time required for the first molecule of kinase c to encounter the nuclear pore) varies in relation to changes in relative component concentrations. For each set of simulations, we observe a distribution of times of signal transduction, demonstrating the inherent noise caused by the limited number of components in the system. Signaling efficiency (as determined by the median signal transduction time) was greatest when the relative concentration of each diffusible component was equivalent (kinase b = kinase c = kinase d = 20). Where the concentration of a particular component decreases, we witness a broadening in the distribution of signal transduction times reflecting the increasing influence of a rate limiting step into the system. We next investigated the effect of co-localizing kinase components to a discrete part of the cell by the inclusion of a ‘constraining’ environment (Methods). Figure 2B and C reveal that co-localizing components within the same environment enhances their ability to transmit a signal. The greatest effect was observed when all three cytosolic components were restricted within the constraining environment, leading to a higher proportion of simulations completing within fewer iterations, compared with simulations performed with no constraining influences. Similar behavior was also observed when only two components were co-localized. Of note, the pathway relationship between the components also appeared to have an effect, with the co-localization of components kinase a and kinase c, increasing the speed of signal transduction relative to those simulations in which kinase a and kinase b (or kinase b and kinase c—data not shown) were co-localized. This effect is presumably related to the fact that both components kinase a and kinase c directly interact with the membrane receptor and nuclear pore, respectively.
|
Finally, we examined the influence of three different architectures of molecular crowding environments. These consisted of filled (inaccessible) lattice sites that were:
- randomly distributed—‘single particles’;
- joined to form rods randomly connecting two sites on the outer membrane—‘transverse rods’ or
- joined to form rods that connected a site on the nuclear membrane to a site on the outer membrane—‘radial rods’ (rods were chosen as a coarse approximation of cytoskeletal fibres and associated proteins).
Results from simulations where 20% of the cell volume was filled with crowding particles (Supplementary Fig. 1). Each crowding environment was found to increase the time of signal transduction. The greatest effect was observed for simulations involving transverse rods (32% of simulations took >10 000 iterations to complete) compared to radial rods (19%), single particles (16%) and no crowding (11%).
These results show that compartmentalization of common pathway components can play a major role in enhancing signal transduction efficiency, suggesting an alternative mechanism to complex formation (Kholodenko, 2003). Furthermore, while the precise relationships with component diffusion and compartmentalization are unclear, our initial studies demonstrate that molecular crowding impacts the efficiency of signal transduction, and that the degree of impact may be related to its organization.
3.2 Case study 2—a metabolic pathway
To optimize their efficiency, many metabolic pathways are highly organized, through sub-cellular compartmentalization and/or complex formation, allowing the channeling of substrates to their target enzymes (Jorgensen et al., 2005). Such metabolic channeling may provide additional advantages including the rapid turnover of potentially toxic intermediates and tighter control over possible metabolic cross talk. Here, we investigate the ability of Cell++ to model a simple metabolic pathway and explore the potential influence of compartmentalization.
Based on Vmax values quoted previously (Teusink et al., 2000), we initially noted that the rate of substrate diffusion effectively negated the effects of compartmentalization (data not shown). However, given the conflict with other studies (Schomburg et al., 2002) (Systems and Methods), simulations were performed in which the value of Vmax was increased by 3–5 orders of magnitude relative to those listed in Supplementary Table 1. Figure 3 shows a view from a typical simulation.
|
Figure 4 presents a series of graphs showing temporal changes in substrate concentrations as a function of enzyme localization and relative reaction speed (increase in Vmax). In the first set of graphs (Fig. 4A), each enzyme is free to diffuse throughout the system. Due to their uniform distribution, these simulations are analogous to ODE models representing homogeneously distributed enzymes and substrates. Consequently, the rate of pyruvate production increases in proportion to the relative increase in Vmax. This is also shown in Supplementary Figure 1, in which T50 values (the time at which 50% of the initial 1,3-bisphosphoglycerate is converted into pyruvate) provide a measure of the overall reaction efficiency. Examining Figure 5A more in detail, the high Keq and Vmax associated with phosphoglycerate kinase results in the rapid production of 3-phosphoglycerate. Subsequent conversion to pyruvate proceeds at a constant, albeit relatively lower, rate (presumably related to the lower Vmax values associated with the intermediate reactions).
|
Constraint of enzymes led to less intuitive patterns of behavior, with a reduction in the rate of pyruvate production dependent upon their pathway relationships (Fig. 4 and Supplementary Fig. 2). For example, whereas the T50 value at a relative speed of reaction of 10 000 was 33 ms for the simulation involving no enzyme constraint, the values for simulations in which: all four enzymes; only phosphoglycerate kinase and pyruvate kinase; and phosphoglycerate kinase and enolase were constrained were 49.5, 99.5 and 198.5 ms, respectively. However, at the lower relative speeds of reaction of 1000, co-localization of all four enzymes resulted in a lower T50 values compared with their free diffusion (296.5 versus 310 ms, respectively).
Closer examination of these data reveals two further relationships, reminiscent of the phenomenon of metabolic channeling. First, the increase in catalytic efficiency, due to co-localization, is greatest for reactions with low equilibrium constants. Second, in the localized absence of the subsequent enzyme in the pathway, intermediates tend to accumulate. The first effect can be seen from comparisons of Figure 4B and E with 4C and D. Since phosphoglycerate mutase has a relatively low equilibrium constant (0.19), its reaction rate is critically dependent on the local concentration of its product, 2-phosphoglycerate. In the simulations presented in 4B and 4E, enolase is constrained to the same locations as phosphoglycerate mutase, ensuring that 2-phosphoglycerate is rapidly converted to the next substrate in the pathway (phosphoenolpyruvate). This removes the kinetic constraint on the first reaction step, resulting in a decrease in the concentration of 3-phosphoglycerate and increase in overall reaction efficiency. In contrast, the simulations presented in Figure 4C and D were performed by constraining either phosphoglycerate mutase or enolase to the central lattice site. Consequently, due to a much reduced concentration of the other enzyme in the same vicinity as the constrained enzyme, the concentration of 3-phosphoglycerate remains relatively high throughout the simulations and the overall reaction speed is not significantly increased by increasing the relative speed of reaction. The second interesting observation (accumulation of intermediates in the localized absence of appropriate enzymes) arises from comparisons of Figure 4C and E with 4B and D. In the former, a substrate (2-phosphoglycerate and phosphoenol pyruvate, respectively) produced within the bulk cytoplasm must diffuse to the central section of the simulation environment in which, enolase and pyruvate kinase are constrained, respectively, before conversion to the next substrate intermediate. This leads to an increase in their relative concentration. In experiments used to produce the latter set of graphs, intermediate substrates are either processed immediately by co-localized enzymes located in the central lattice site (Fig. 4B), or diffuse a short distance to the bulk cytoplasm where they are rapidly converted to pyruvate (Fig. 4D). Finally, it is worth mentioning that a lower initial concentration of 1,3-bisphosphoglycerate (1 µM total concentration) resulted in qualitatively similar behavior to previous simulations (Fig. 4F and G).
These results demonstrate that co-localization can increase pathway efficiency, the effect being related to both the pathway relationships of the co-localized enzymes and their kinetic parameters (equilibrium constants). Furthermore, enzyme co-localization also led to a decrease in the accumulation of intermediate substrates. These findings are consistent with the phenomenon of metabolic channeling (Jorgensen et al., 2005).
3.3 Case study 3—intracellular calcium signaling
Ionized calcium (Ca2+) is an important second messenger molecule involved in a number of diverse cellular processes including cell differentiation, growth and metabolism (Brini, 2003; Webb and Miller, 2003). Propagation of intracellular Ca2+ signals is mediated by a Ca2+-induced Ca2+-release mechanism (CICR), resulting in the formation of a so called Ca2+ wave (Berridge et al., 2003). For example, in muscle cells, calcium-release channels populate the surface of the sarcoplasmic reticulum (SR). Upon a localized increase in Ca2+ concentration they open, releasing more stored Ca2+ from intraluminal stores. Diffusion of this released Ca2+ activates more release channels, resulting in a wave of Ca2+ propagating over the surface of the SR. The properties of these waves (e.g. speed and amplitude) are thought to be modulated by the kinetic properties of uptake and release channels (ATP-dependent Ca2+ pumps; and inositol triphosphate and ryanodine receptors, respectively) as well as their spatial distribution (Vermassen et al., 2003; Podhaisky and Wussling, 2004). Due to limitations in optical resolution, the precise relationships between spatial organization and CICR remain unclear. Here, we applied Cell++ to assess the impact of three types of release channel organization on a model of CICR (see Supplementary Fig. 3 for a screen shot of a typical simulation).
Results from these experiments qualitatively demonstrate how release channel receptor concentration and spatial organization can affect Ca2+ wave propagation (Fig. 5). For example, for each distributions of channels, there is an exponential relationship between the speed of the wave and the number of Ca2+ channels. In terms of spatial organization, simulations involving the random placement of channels produced similar results to those in which channels were uniformly distributed. However, clustering of channels reduced both the probability that a wave successfully propagated and (with the exception of simulations involving only 10 channels) its average speed. These effects are presumably a consequence of the larger gaps in the CICR network generated by the clustering of channels. It is interesting to note, that for those simulations involving only 10 channels for which wave propagation was successful, the clustering of receptors actually led to faster wave propagation (4476 iterations versus 4872 and 4787 iterations for the uniform and random distributions, respectively). Together these results show that the spatial distribution of ryanodine receptors affects both the speed and successful propagation of CICR. It remains an interesting conjecture to propose that their organization may at least in part, be responsible for generating different properties of Ca2+ waves (e.g. relative speeds, amplitude and frequency). By incorporating additional factors such as Ca2+ uptake (e.g. via SERCA pumps), we suggest that Cell++ is ideally suited to explore these relationships in more detail.
|
|
4 DISCUSSION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSIONREFERENCES |
|---|
With the availability of increasingly large datasets detailing cellular components and their interactions, there is a growing need to develop novel tools capable of integrating these data to provide intuitive dynamic models of biochemical pathways. Furthermore, given the potential impact of spatial organization (Vermassen et al., 2003; Ovadi and Saks, 2004; Jorgensen et al., 2005), it is imperative that these tools allow the incorporation of such factors. Cell++ has been designed to meet these needs. In a unique approach, Cell++ combines Brownian dynamics to model ‘large’ discrete molecules, such as enzymes, with a cellular automata engine to describe small molecule behavior. By adopting this dual approach, simulations involving large numbers of molecules can be performed without significantly compromising the computational overhead. As an example, the metabolic pathway simulations described here, which involved 4000 discrete enzymes and
1.25 x 10–19 mols (or 92 000 molecules) of substrate, were typically completed within 2–4 h by a single Intel (2.4 GHz) processor. Cell++ thus provides the user with a fast and flexible simulation environment, which together with its real time visualization interfaces, is expected to find wide application within both research and educational contexts. Below are highlighted some key features:- Flexibility. Cell++ provides a flexible simulation environment that models both individual discrete molecules and larger populations of small molecule concentrations.
- A visualization interface. The real time 3D visualization of simulations by Cell++, permits the user to readily investigate the behavior of the system and gain insights that are not readily obtained via numerical outputs.
- A spatial context. One of our main goals in creating Cell++, was to provide a tool that could readily investigate the influence of spatial factors on cellular processes. By adopting a 3D lattice, it is trivial to define different spatial environments that either restrict component movement or alter their behavior.
- Speed. We noted that for simulations (e.g. signaling pathways), which may be prone to stochastic influences, it is important to examine the output of many similar simulations from a statistical perspective. As observed from our simulations of signaling and metabolic pathways, Cell++ can provide results from several thousand simulations in only a few hours on a single Linux workstation.
- Accessibility. Cell++ is freely available under the GNU general public license (GPL) agreement (http://www.gnu.org/copyleft/gpl.html). A beta version of the software together with an overview of the software, some applications and a tutorial is available at http://www.compsysbio.org/CellSim/. The simulations involving calcium waves were performed using a slightly modified version of the available source code. Details of its implementation are available upon request to the authors.
Previously, a number of biochemical simulation environments have been developed and applied to the study of particular classes of systems, from simple signaling pathways to more complex systems, such as the release of neurotransmitters in synapses and membrane growth (Tomita et al., 1999; Kierzek, 2002; Ander et al., 2004; Broderick et al., 2004; Hattne et al., 2005; Coggan et al., 2005). Among the better known of these platforms is the ODE-based package ‘E-CELL’ (Tomita et al., 1999). As a way of comparison between E-CELL and Cell++, we applied the former to perform an analogous simulation involving the glycolysis metabolic pathway presented above, in which enzymes were either free to roam throughout the simulation environment or constrained to a smaller compartment (Supplementary Methods for further details of the E-CELL implementation).
Qualitatively the overall results from these simulations, in terms of changes in metabolite concentrations, displayed similar behaviors to those obtained with Cell++ (Supplementary Figs 4 and 5). For example, results for simulations in which enzymes were free to roam the entire system, displayed very similar behaviors (both qualitative and quantitative) for the two platforms (Supplementary Fig. 4A). This is perhaps unsurprising as due to the large number of enzymes used; this system represents the well-stirred reaction vessel model implicit with E-CELL. However, the incorporation of spatial constraints led to differences in the shapes and sizes of the resulting graphs of metabolite concentrations (Supplementary Fig. 4B and C). This reflects differences in the way in which intercompartmental diffusion is handled by the two platforms and we found that by altering the diffusion constant in the E-CELL simulations (Supplementary Fig. 4C), we were able to dramatically alter metabolite turnover. We also examined the effect of altering Vmax and noted that qualitatively, E-CELL and Cell++ again resulted in similar behaviors (Supplementary Fig. 5).
Given the large number of enzymes and metabolites used in these studies, it is therefore reassuring to note that Cell++ simulations were consistent with those derived using E-CELL. We suggest therefore that the since the E-CELL environment provides an intuitive interface for construction of models of well mixed systems, that E-CELL remains an ideal tool for the study of such systems. However, as we noted in setting up these systems initially, it may be difficult to accurately implement diffusion parameters for complex processes involving more detailed spatial architectures. Hence the explicit diffusion model employed by Cell++ makes it more suitable for the study of spatially complex systems. Furthermore, it may also be preferable to employ Cell++ for the study of more stochastic systems composed of e.g. limited numbers of enzymes and/or substrates (e.g. simulations of calcium waves).
While a comprehensive survey of each of the systems presented here is clearly outside the current scope of this paper, findings from our initial studies nonetheless highlight the potential importance of spatial organization and its impact on the evolution and engineering of biochemical pathways.
Given the vast array of molecular processes operating within a cell, the number of proteins assigned to any particular pathway is clearly limited. Previous studies of protein expression in yeast found that the number of each type of protein in a cell can vary from 50 to several million (Ghaemmaghami et al., 2003). While protein expression may be highly sensitive to e.g. changes in the environment and the cell cycle, it is presumed that the relative expression of pathway components has been optimized to ensure their efficient operation. However, the spatial organization of pathways provides an additional mechanism to increase pathway efficiency and thereby limit the allocation of protein resources associated with maintaining a diverse set of biochemical pathways. In addition, component co-localization potentially offers the additional benefit of minimizing pathway cross talk (Ahmed et al., 2005). For example, it is well known that many signaling pathways share common components [e.g. the gene products of STE7 and STE11 are used in at least three distinct MAPK based signaling pathways in yeast (Flautner et al., 2005)]. By compartmentalizing overlapping pathways, signal fidelities and specificities are maintained (Kondoh et al., 2005). This leads to the intriguing conjecture that the evolution of spatial organization, as a means of increasing pathway efficiency, allows the creation of novel pathways by the adoption and spatial segregation of components from existing pathways. Such a mechanism has the potential to increase the repertoire of signaling systems without compromising the need for a significant investment in additional resources. In metabolic pathways, compartmentalization of common enzymes can lead to the metabolic channeling of substrates through particular pathways within the metabolic network (Jorgensen et al., 2005). This provides a degree of control over metabolic cross talk that can arise from the sharing of substrate intermediates between different pathways. Thus, in addition to changes in protein expression, changes in the co-localization of enzymes may also underlie a cells ability rapidly to switch pathways, perhaps in response to changes in its local environment (Heilmann et al., 2004). It is expected that more detailed studies of larger systems of interacting metabolic pathways will help uncover additional relationships between pathway architectures, reaction kinetics and spatial organization, thereby increasing the potential of targeted metabolic engineering.
Future directions of Cell++ include the development of additional documentation and improvement of the interfaces used to configure simulations. We are currently investigating integration of the systems biology markup language—SBML—into Cell++ (Hucka et al., 2003). In addition to simplifying the process of system configuration through the use of existing SBML model editors such as CellDesigner (Funahasi et al., 2003), this would offer the additional benefit of providing a mechanism to readily interconvert models between Cell++ and other simulation tools.
|
Acknowledgments |
|---|
The authors would like to thank the Center for Computational Biology, Hospital for Sick Children, for computing facilities and Alan Davidson and Shoshana Wodak for comments during manuscript preparation. C.S. is supported by the Hospital for Sick Children (Toronto, ON, Canada) Research Training Centre. J.P. is supported by a Canadian Institute of Health Research New Investigators Award. This work was funded by a discovery grant from the Natural Sciences and Engineering Research Council of Canada.
Conflict of Interest: none declared.
|
FOOTNOTES |
|---|
The authors wish it to be known that, in their opinion, the first two authors should be regarded as joint First Authors. 
Associate Editor: Martin Bishop
Received on June 1, 2006; revised on September 21, 2006; accepted on September 27, 2006
|
REFERENCES |
|---|
|
TOPABSTRACT1 INTRODUCTION2 METHODS3 RESULTS4 DISCUSSIONREFERENCES |
|---|
Ahmed, Z., et al. (2005) Distinct spatial and temporal distribution of ZAP70 and Lck following stimulation of interferon and T-cell receptors. J. Mol. Biol, . 353, 1001–1010[CrossRef][Web of Science][Medline].
Ander, M., et al. (2004) SmartCell, a framework to simulate cellular processes that combines stochastic approximation with diffusion and localization: analysis of simple networks. Syst. Biol, . 1, 129–138[CrossRef].
Berridge, M.J., et al. (2003) Calcium signaling: dynamics, homeostasis and remodeling. Nat. Rev. Mol. Cell. Biol, . 4, 517–529[CrossRef][Web of Science][Medline].
Brini, M. (2003) Ca(2+) signaling in mitochondria: mechanism and role in physiology and pathology. Cell. Calcium, 34, 399–405[CrossRef][Web of Science][Medline].
Broderick, G., et al. (2004) A life-like virtual cell membrane using discrete automata. In Silico Biol, . 5, 0016.
Coggan, J.S., et al. (2005) Evidence for ectopic neurotransmission at a neuronal synapse. Science, 309, 446–451
Dayel, M.J., et al. (1999) Diffusion of green fluorescent protein in the aqueous-phase lumen of endoplasmic reticulum. Biophys J, . 76, 2843–2851.
Elowitz, M., et al. (1999) Protein mobility in the cytoplasm of the Escherichia coli. J. Bacteriol, . 181, 197–203
Flatauer, L.J., et al. (2005) Mitogen-activated protein kinases with distinct requirements for Ste5 scaffolding influence signaling specificity in Saccharomyces cerevisiae. Mol. Cell. Biol, . 25, 1793–1803
Funahashi, A., et al. (2003) CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. Biosilico, 1, 159–162[CrossRef].
Ghaemmaghami, S., et al. (2003) Global analysis of protein expression in yeast. Nature, 425, 737–741[CrossRef][Medline].
Hattne, J., et al. (2005) Stochastic reaction-diffusion simulation with MesoRD. Bioinformatics, 21, 2923–2924
Heilmann, I., et al. (2004) Switching desaturase enzyme specificity by alternate subcellular targeting. Proc. Natl Acad. Sci. USA, 101, 10266–10271
Hermjakob, H., et al. (2004) IntAct—an open source molecular interaction database. Nucleic Acids Res, . 32, D452–D455
Hucka, M., et al. (2003) The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics, 19, 524–531
Ingber, D.E., et al. (1994) Cellular tensegrity: exploring how mechanical changes in the cytoskeleton regulate cell growth, migration, and tissue pattern during morphogenesis. Int. Rev. Cytol, . 150, 173–224[Web of Science][Medline].
Ingber, D.E. (2003) Tensegrity I. Cell structure and hierarchical systems biology. J. Cell. Sci, . 116, 1157–1173
Jorgensen, K., et al. (2005) Metabolon formation and metabolic channeling in the biosynthesis of plant natural products. Curr. Opin. Plant Biol, . 8, 280–291[CrossRef][Web of Science][Medline].
Kholodenko, B.N. (2003) Four-dimensional organization of protein kinase signaling cascades: the roles of diffusion, endocytosis and molecular motors. J. Exp. Biol, . 206, 2073–2082
Kierzek, A.M. (2002) STOCKS: STOChastic Kinetic Simulations of biochemical systems with Gillespie algorithm. Bioinformatics, 18, 470–481
Kondoh, K., Torii, S., Nishida, E. (2005) Control of MAP kinase signaling to the nucleus. Chromosoma, 114, 86–91[CrossRef][Web of Science][Medline].
Monastyrskaya, K., et al. (2005) The NK1 receptor localizes to the plasma membrane microdomains, and its activation is dependent on lipid raft integrity. J. Biol. Chem, . 280, 7135–7146
Ovadi, J. and Saks, V. (2004) On the origin of intracellular compartmentation and organized metabolic systems. Mol. Cell. Biochem, . 5–12.
Podhaisky, H. and Wussling, M.H. (2004) The velocity of calcium waves is expected to depend non-monotoneously on the density of the calcium release units. Mol. Cell. Biochem, . 387–390.
Salwinski, L., et al. (2004) The database of interacting proteins: 2004 update. Nucleic Acids Res, . 32, D449–D451
Schomburg, I., et al. (2002) BRENDA, enzyme data and metabolic information. Nucleic Acids Res, . 30, 47–49
Takahashi, K., et al. (2005) Space in systems biology of signaling pathways–towards intracellular molecular crowding in silico. FEBS Lett, . 579, 1783–1788[CrossRef][Web of Science][Medline].
Teusink, B., et al. (2000) Can yeast glycolysis be understood in terms of in vitro kinetics of the constituent enzymes? Testing biochemistry. Eur. J. Biochem, . 267, 5313–5329[Web of Science][Medline].
Tomita, M., et al. (1999) E-CELL: software environment for whole-cell simulation. Bioinformatics, 15, 72–84
Verkman, A.S. (2002) Solute and macromolecule diffusion in cellular aqueous compartments. Trends Biochem. Sci, . 27, 27–33[CrossRef][Web of Science][Medline].
Vermassen, E., et al. (2003) Microtubule-dependent redistribution of the type-1 inositol 1,4,5-triphosphate receptor in A7r5 smooth muscle cells. J. Cell. Sci, . 116, 1269–1277
Webb, S.E. and Miller, A.L. (2003) Calcium signaling during embryonic development. Nat. Rev. Mol. Cell. Biol, . 4, 539–551[CrossRef][Web of Science][Medline].
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||