Structure-based kinetic models of modular signaling protein function: focus on Shp2

An integrated mechanistic and data-driven computational model predicts cell responses to high- and low-affinity EGFR ligands

The biophysical properties of ligand binding heavily influence the ability of receptors to specify cell fates. Understanding the rules by which ligand binding kinetics impact cell phenotype is challenging, however, because of the coupled information transfers that occur from receptors to downstream signaling effectors and from effectors to phenotypes. Here, we address that issue by developing an integrated mechanistic and data-driven computational modeling platform to predict cell responses to different ligands for the epidermal growth factor receptor (EGFR). Experimental data for model training and validation were generated using MCF7 human breast cancer cells treated with the high- and low-affinity ligands epidermal growth factor (EGF) and epiregulin (EREG), respectively. The integrated model captures the unintuitive, concentration-dependent abilities of EGF and EREG to drive signals and phenotypes differently, even at similar levels of receptor occupancy. For example, the model correctly predicts the dominance of EREG over EGF in driving a cell differentiation phenotype through AKT signaling at intermediate and saturating ligand concentrations and the ability of EGF and EREG to drive a broadly concentration-sensitive migration phenotype through cooperative ERK and AKT signaling. Parameter sensitivity analysis identifies EGFR endocytosis, which is differentially regulated by EGF and EREG, as one of the most important determinants of the alternative phenotypes driven by different ligands. The integrated model provides a new platform to predict how phenotypes are controlled by the earliest biophysical rate processes in signal transduction and may eventually be leveraged to understand receptor signaling system performance depends on cell context.

One-sentence summary

Integrated kinetic and data-driven EGFR signaling model identifies the specific signaling mechanisms that dictate cell responses to EGFR activation by different ligands.

Dissecting protein tyrosine phosphatase signaling by engineered chemogenetic control of its activity

Protein tyrosine phosphatases (PTPases) are critical mediators of dynamic cell signaling. A tool capable of identifying transient signaling events downstream of PTPases is essential to understand phosphatase function on a physiological time scale. We report a broadly applicable protein engineering method for allosteric regulation of PTPases. This method enables dissection of transient events and reconstruction of individual signaling pathways. Implementation of this approach for Shp2 phosphatase revealed parallel MAPK and ROCK II dependent pathways downstream of Shp2, mediating transient cell spreading and migration. Furthermore, we show that the N-SH2 domain of Shp2 regulates MAPK-independent, ROCK II-dependent cell migration. Engineered targeting of Shp2 activity to different protein complexes revealed that Shp2-FAK signaling induces cell spreading whereas Shp2-Gab1 or Shp2-Gab2 mediates cell migration. We identified specific transient morphodynamic processes induced by Shp2 and determined the role of individual signaling pathways downstream of Shp2 in regulating these events. Broad application of this approach is demonstrated by regulating PTP1B and PTP-PEST phosphatases.

Understanding the Dynamics of T-Cell Activation in Health and Disease Through the Lens of Computational Modeling

T cells in the immune system are activated by binding to foreign peptides (from an external pathogen) or mutant peptide (derived from endogenous proteins) displayed on the surface of a diseased cell. This triggers a series of intracellular signaling pathways, which ultimately dictate the response of the T cell. The insights from computational models have greatly improved our understanding of the mechanisms that control T-cell activation. In this review, we focus on the use of ordinary differential equation–based mechanistic models to study T-cell activation. We highlight several examples that demonstrate the models’ utility in answering specific questions related to T-cell activation signaling, from antigen discrimination to the feedback mechanisms that initiate transcription factor activation. In addition, we describe other modeling approaches that can be combined with mechanistic models to bridge time scales and better understand how intracellular signaling events, which occur on the order of seconds to minutes, influence phenotypic responses of T-cell activation, which occur on the order of hours to days. Overall, through concrete examples, we emphasize how computational modeling can be used to enable the rational design and optimization of immunotherapies.

Analysis of Fn14-NF-κB signaling response dynamics using a mechanistic model

Fn14 is a transmembrane receptor protein belonging to the tumor necrosis factor receptor (TNFR) superfamily. Many experimental reports have shown that crosslinking of the receptor by its extracellular ligand TWEAK induces prolonged activation of transcription factor NF-κB. This behavior is distinct from TNF-α receptor, which is a more well-characterized member of the TNFR family. TNF-α receptor, despite sharing many similar molecular interactions with Fn14, only transiently activates NF-κB in response to TNF-α stimulation. Here, we investigate molecular mechanisms that enable Fn14 to display such distinctive behavior. In particular, we focus on two specific features of the Fn14 pathway that potentially give rise to a positive feedback regulation and differentiate it from the TNF-α receptor signaling. By developing a mechanistic model, we analyze how these features may determine the dynamics of an Fn14-NF-κB response. Our analysis reveals that stimulation of Fn14 by TWEAK may generate highly non-linear dynamics, including stable limit cycles and bistable responses. The type of response depends both on the strength and duration of a TWEAK signal. Our predictions and analyses also show that the molecular interactions underlying the positive feedback explain the prolonged activation of NF-κB under certain parameter regimes. In light of the model predictions, we propose possible deregulations of Fn14 leading to its overexpression in solid tumors and tissue injuries.

Mathematical modeling identifies Lck as a potential mediator for PD-1 induced inhibition of early TCR signaling

Programmed cell death-1 (PD-1) is an inhibitory immune checkpoint receptor that negatively regulates the functioning of T cell. Although the direct targets of PD-1 were not identified, its inhibitory action on the TCR signaling pathway was known much earlier. Recent experiments suggest that the PD-1 inhibits the TCR and CD28 signaling pathways at a very early stage ─ at the level of phosphorylation of the cytoplasmic domain of TCR and CD28 receptors. Here, we develop a mathematical model to investigate the influence of inhibitory effect of PD-1 on the activation of early TCR and CD28 signaling molecules. Proposed model recaptures several quantitative experimental observations of PD-1 mediated inhibition. Model simulations show that PD-1 imposes a net inhibitory effect on the Lck kinase. Further, the inhibitory effect of PD-1 on the activation of TCR signaling molecules such as Zap70 and SLP76 is significantly enhanced by the PD-1 mediated inhibition of Lck. These results suggest a critical role for Lck as a mediator for PD-1 induced inhibition of TCR signaling network. Multi parametric sensitivity analysis explores the effect of parameter uncertainty on model simulations.

Evaluation of Parallel Tempering to Accelerate Bayesian Parameter Estimation in Systems Biology

Models of biological systems often have many unknown parameters that must be determined in order for model behavior to match experimental observations. Commonly-used methods for parameter estimation that return point estimates of the best-fit parameters are insufficient when models are high dimensional and under-constrained. As a result, Bayesian methods, which treat model parameters as random variables and attempt to estimate their probability distributions given data, have become popular in systems biology. Bayesian parameter estimation often relies on Markov Chain Monte Carlo (MCMC) methods to sample model parameter distributions, but the slow convergence of MCMC sampling can be a major bottleneck. One approach to improving performance is parallel tempering (PT), a physics-based method that uses swapping between multiple Markov chains run in parallel at different temperatures to accelerate sampling. The temperature of a Markov chain determines the probability of accepting an unfavorable move, so swapping with higher temperatures chains enables the sampling chain to escape from local minima. In this work we compared the MCMC performance of PT and the commonly-used Metropolis-Hastings (MH) algorithm on six biological models of varying complexity. We found that for simpler models PT accelerated convergence and sampling, and that for more complex models, PT often converged in cases MH became trapped in non-optimal local minima. We also developed a freely-available MATLAB package for Bayesian parameter estimation called PTEMPEST (http://github.com/RuleWorld/ptempest), which is closely integrated with the popular BioNetGen software for rule-based modeling of biological systems.

Topological and kinetic determinants of the modal matrices of dynamic models of metabolism

Large-scale kinetic models of metabolism are becoming increasingly comprehensive and accurate. A key challenge is to understand the biochemical basis of the dynamic properties of these models. Linear analysis methods are well-established as useful tools for characterizing the dynamic response of metabolic networks. Central to linear analysis methods are two key matrices: the Jacobian matrix (J) and the modal matrix (M-1) arising from its eigendecomposition. The modal matrix M-1 contains dynamically independent motions of the kinetic model near a reference state, and it is sparse in practice for metabolic networks. However, connecting the structure of M-1 to the kinetic properties of the underlying reactions is non-trivial. In this study, we analyze the relationship between J, M-1, and the kinetic properties of the underlying network for kinetic models of metabolism. Specifically, we describe the origin of mode sparsity structure based on features of the network stoichiometric matrix S and the reaction kinetic gradient matrix G. First, we show that due to the scaling of kinetic parameters in real networks, diagonal dominance occurs in a substantial fraction of the rows of J, resulting in simple modal structures with clear biological interpretations. Then, we show that more complicated modes originate from topologically-connected reactions that have similar reaction elasticities in G. These elasticities represent dynamic equilibrium balances within reactions and are key determinants of modal structure. The work presented should prove useful towards obtaining an understanding of the dynamics of kinetic models of metabolism, which are rooted in the network structure and the kinetic properties of reactions.

A multiscale modeling study of particle size effects on the tissue penetration efficacy of drug-delivery nanoparticles

BACKGROUND

Particle size is a key parameter for drug-delivery nanoparticle design. It is believed that the size of a nanoparticle may have important effects on its ability to overcome the transport barriers in biological tissues. Nonetheless, such effects remain poorly understood. Using a multiscale model, this work investigates particle size effects on the tissue distribution and penetration efficacy of drug-delivery nanoparticles.

RESULTS

We have developed a multiscale spatiotemporal model of nanoparticle transport in biological tissues. The model implements a time-adaptive Brownian Dynamics algorithm that links microscale particle-cell interactions and adhesion dynamics to tissue-scale particle dispersion and penetration. The model accounts for the advection, diffusion, and cellular uptakes of particles. Using the model, we have analyzed how particle size affects the intra-tissue dispersion and penetration of drug delivery nanoparticles. We focused on two published experimental works that investigated particle size effects in in vitro and in vivo tissue conditions. By analyzing experimental data reported in these two studies, we show that particle size effects may appear pronounced in an in vitro cell-free tissue system, such as collagen matrix. In an in vivo tissue system, the effects of particle size could be relatively modest. We provide a detailed analysis on how particle-cell interactions may determine distribution and penetration of nanoparticles in a biological tissue.

CONCLUSION

Our work suggests that the size of a nanoparticle may play a less significant role in its ability to overcome the intra-tissue transport barriers. We show that experiments involving cell-free tissue systems may yield misleading observations of particle size effects due to the absence of advective transport and particle-cell interactions.

Biophysical assay for tethered signaling reactions reveals tether-controlled activity for the phosphatase SHP-1

Tethered enzymatic reactions are ubiquitous in signaling networks but are poorly understood. A previously unreported mathematical analysis is established for tethered signaling reactions in surface plasmon resonance (SPR). Applying the method to the phosphatase SHP-1 interacting with a phosphorylated tether corresponding to an immune receptor cytoplasmic tail provides five biophysical/biochemical constants from a single SPR experiment: two binding rates, two catalytic rates, and a reach parameter. Tether binding increases the activity of SHP-1 by 900-fold through a binding-induced allosteric activation (20-fold) and a more significant increase in local substrate concentration (45-fold). The reach parameter indicates that this local substrate concentration is exquisitely sensitive to receptor clustering. We further show that truncation of the tether leads not only to a lower reach but also to lower binding and catalysis. This work establishes a new framework for studying tethered signaling processes and highlights the tether as a control parameter in clustered receptor signaling.

Shp2 Associates with and Enhances Nephrin Tyrosine Phosphorylation and Is Necessary for Foot Process Spreading in Mouse Models of Podocyte Injury

In most forms of glomerular diseases, loss of size selectivity by the kidney filtration barrier is associated with changes in the morphology of podocytes. The kidney filtration barrier is comprised of the endothelial lining, the glomerular basement membrane, and the podocyte intercellular junction, or slit diaphragm. The cell adhesion proteins nephrin and neph1 localize to the slit diaphragm and transduce signals in a Src family kinase Fyn-mediated tyrosine phosphorylation-dependent manner. Studies in cell culture suggest nephrin phosphorylation-dependent signaling events are primarily involved in regulation of actin dynamics and lamellipodium formation. Nephrin phosphorylation is a proximal event that occurs both during development and following podocyte injury. We hypothesized that abrogation of nephrin phosphorylation following injury would prevent nephrin-dependent actin remodeling and foot process morphological changes. Utilizing a biased screening approach, we found nonreceptor Src homology 2 (sh2) domain-containing phosphatase Shp2 to be associated with phosphorylated nephrin. We observed an increase in nephrin tyrosine phosphorylation in the presence of Shp2 in cell culture studies. In the human glomerulopathies minimal-change nephrosis and membranous nephropathy, there is an increase in Shp2 phosphorylation, a marker of increased Shp2 activity. Mouse podocytes lacking Shp2 do not develop foot process spreading when subjected to podocyte injury in vivo using protamine sulfate or nephrotoxic serum (NTS). In the NTS model, we observed a lack of foot process spreading in mouse podocytes with Shp2 deleted and smaller amounts of proteinuria. Taken together, these results suggest that Shp2-dependent signaling events are necessary for changes in foot process structure and function following injury.

Phosphatase specificity and pathway insulation in signaling networks

Phosphatases play an important role in cellular signaling networks by regulating the phosphorylation state of proteins. Phosphatases are classically considered to be promiscuous, acting on tens to hundreds of different substrates. We recently demonstrated that a shared phosphatase can couple the responses of two proteins to incoming signals, even if those two substrates are from otherwise isolated areas of the network. This finding raises a potential paradox: if phosphatases are indeed highly promiscuous, how do cells insulate themselves against unwanted crosstalk? Here, we use mathematical models to explore three possible insulation mechanisms. One approach involves evolving phosphatase KM values that are large enough to prevent saturation by the phosphatase's substrates. Although this is an effective method for generating isolation, the phosphatase becomes a highly inefficient enzyme, which prevents the system from achieving switch-like responses and can result in slow response kinetics. We also explore the idea that substrate degradation can serve as an effective phosphatase. Assuming that degradation is unsaturatable, this mechanism could insulate substrates from crosstalk, but it would also preclude ultrasensitive responses and would require very high substrate turnover to achieve rapid dephosphorylation kinetics. Finally, we show that adaptor subunits, such as those found on phosphatases like PP2A, can provide effective insulation against phosphatase crosstalk, but only if their binding to substrates is uncoupled from their binding to the catalytic core. Analysis of the interaction network of PP2A's adaptor domains reveals that although its adaptors may isolate subsets of targets from one another, there is still a strong potential for phosphatase crosstalk within those subsets. Understanding how phosphatase crosstalk and the insulation mechanisms described here impact the function and evolution of signaling networks represents a major challenge for experimental and computational systems biology.

Analytical reduction of combinatorial complexity arising from multiple protein modification sites

Combinatorial complexity is a major obstacle to ordinary differential equation (ODE) modelling of biochemical networks. For example, a protein with 10 sites that can each be unphosphorylated, phosphorylated or bound to adaptor protein requires 3(10) ODEs. This problem is often dealt with by making ad hoc assumptions which have unclear validity and disallow modelling of site-specific dynamics. Such site-specific dynamics, however, are important in many biological systems. We show here that for a common biological situation where adaptors bind modified sites, binding is slow relative to modification/demodification, and binding to one modified site hinders binding to other sites, for a protein with n modification sites and m adaptor proteins the number of ODEs needed to simulate the site-specific dynamics of biologically relevant, lumped bound adaptor states is independent of the number of modification sites and equal to m + 1, giving a significant reduction in system size. These considerations can be relaxed considerably while retaining reasonably accurate descriptions of the true system dynamics. We apply the theory to model, using only 11 ODEs, the dynamics of ligand-induced phosphorylation of nine tyrosines on epidermal growth factor receptor (EGFR) and primary recruitment of six signalling proteins (Grb2, PI3K, PLCγ1, SHP2, RasA1 and Shc1). The model quantitatively accounts for experimentally determined site-specific phosphorylation and dephosphorylation rates, differential affinities of binding proteins for the phosphorylated sites and binding protein expression levels. Analysis suggests that local concentration of site-specific phosphatases such as SHP2 in membrane subdomains by a factor of approximately 10(7) is critical for effective site-specific regulation. We further show how our framework can be extended with minimal effort to consider binding cooperativity between Grb2 and c-Cbl, which is important for receptor trafficking. Our theory has potentially broad application to reduce combinatorial complexity and allow practical simulation of a variety ODE models relevant to systems biology and pharmacology applications to allow exploration of key aspects of complexity that control signal flux.

Systems biology of IL-6, IL-12 family cytokines

Interleukin-6-type cytokines play important roles in the communication between cells of multicellular organisms. They are involved in the regulation of complex cellular processes such as proliferation and differentiation and act as key player during inflammation and immune response. A major challenge is to understand how these complex non-linear processes are connected and regulated. Systems biology approaches are used to tackle this challenge in an iterative process of quantitative experimental and mathematical analyses. Here we review quantitative experimental studies and systems biology approaches dealing with the function of Interleukin-6-type cytokines in physiological and pathophysiological conditions. These approaches cover the analyses of signal transduction on a cellular level up to pharmacokinetic and pharmacodynamic studies on a whole organism level.

Modeling biomolecular site dynamics in immunoreceptor signaling systems

The immune system plays a central role in human health. The activities of immune cells, whether defending an organism from disease or triggering a pathological condition such as autoimmunity, are driven by the molecular machinery of cellular signaling systems. Decades of experimentation have elucidated many of the biomolecules and interactions involved in immune signaling and regulation, and recently developed technologies have led to new types of quantitative, systems-level data. To integrate such information and develop nontrivial insights into the immune system, computational modeling is needed, and it is essential for modeling methods to keep pace with experimental advances. In this chapter, we focus on the dynamic, site-specific, and context-dependent nature of interactions in immunoreceptor signaling (i.e., the biomolecular site dynamics of immunoreceptor signaling), the challenges associated with capturing these details in computational models, and how these challenges have been met through use of rule-based modeling approaches.

Receptor tyrosine kinase (RTK) signalling in the control of neural stem and progenitor cell (NSPC) development

Important developmental responses are elicited in neural stem and progenitor cells (NSPC) by activation of the receptor tyrosine kinases (RTK), including the fibroblast growth factor receptors, epidermal growth factor receptor, platelet-derived growth factor receptors and insulin-like growth factor receptor (IGF1R). Signalling through these RTK is necessary and sufficient for driving a number of developmental processes in the central nervous system. Within each of the four RTK families discussed here, receptors are activated by sets of ligands that do not cross-activate receptors of the other three families, and therefore, their activation can be independently regulated by ligand availability. These RTK pathways converge on a conserved core of signalling molecules, but differences between the receptors in utilisation of signalling molecules and molecular adaptors for intracellular signal propagation become increasingly apparent. Intracellular inhibitors of RTK signalling are widely involved in the regulation of developmental signalling in NSPC and often determine developmental outcomes of RTK activation. In addition, cellular responses of NSPC to the activation of a given RTK may be significantly modulated by signal strength. Cellular propensity to respond also plays a role in developmental outcomes of RTK signalling. In combination, these mechanisms regulate the balance between NSPC maintenance and differentiation during development and in adulthood. Attribution of particular developmental responses of NSPC to specific pathways of RTK signalling becomes increasingly elucidated. Co-activation of several RTK in developing NSPC is common, and analysis of co-operation between their signalling pathways may advance knowledge of RTK role in NSPC development.

Modeling formalisms in Systems Biology

Systems Biology has taken advantage of computational tools and high-throughput experimental data to model several biological processes. These include signaling, gene regulatory, and metabolic networks. However, most of these models are specific to each kind of network. Their interconnection demands a whole-cell modeling framework for a complete understanding of cellular systems. We describe the features required by an integrated framework for modeling, analyzing and simulating biological processes, and review several modeling formalisms that have been used in Systems Biology including Boolean networks, Bayesian networks, Petri nets, process algebras, constraint-based models, differential equations, rule-based models, interacting state machines, cellular automata, and agent-based models. We compare the features provided by different formalisms, and discuss recent approaches in the integration of these formalisms, as well as possible directions for the future.

Guidelines for visualizing and annotating rule-based models

Rule-based modeling provides a means to represent cell signaling systems in a way that captures site-specific details of molecular interactions. For rule-based models to be more widely understood and (re)used, conventions for model visualization and annotation are needed. We have developed the concepts of an extended contact map and a model guide for illustrating and annotating rule-based models. An extended contact map represents the scope of a model by providing an illustration of each molecule, molecular component, direct physical interaction, post-translational modification, and enzyme-substrate relationship considered in a model. A map can also illustrate allosteric effects, structural relationships among molecular components, and compartmental locations of molecules. A model guide associates elements of a contact map with annotation and elements of an underlying model, which may be fully or partially specified. A guide can also serve to document the biological knowledge upon which a model is based. We provide examples of a map and guide for a published rule-based model that characterizes early events in IgE receptor (FcεRI) signaling. We also provide examples of how to visualize a variety of processes that are common in cell signaling systems but not considered in the example model, such as ubiquitination. An extended contact map and an associated guide can document knowledge of a cell signaling system in a form that is visual as well as executable. As a tool for model annotation, a map and guide can communicate the content of a model clearly and with precision, even for large models.

Exact solutions to a spatially extended model of kinase-receptor interaction

B and Mast cells are activated by the aggregation of the immune receptors. Motivated by this phenomena we consider a simple spatially extended model of mutual interaction of kinases and membrane receptors. It is assumed that kinase activates membrane receptors and in turn the kinase molecules bound to the active receptors are activated by transphosphorylation. Such a type of interaction implies positive feedback and may lead to bistability. In this study we apply the Steklov eigenproblem theory to analyze the linearized model and find exact solutions in the case of non-uniformly distributed membrane receptors. This approach allows us to determine the critical value of receptor dephosphorylation rate at which cell activation (by arbitrary small perturbation of the inactive state) is possible. We found that cell sensitivity grows with decreasing kinase diffusion and increasing anisotropy of the receptor distribution. Moreover, these two effects are cooperating. We showed that the cell activity can be abruptly triggered by the formation of the receptor aggregate. Since the considered activation mechanism is not based on receptor crosslinking by polyvalent antigens, the proposed model can also explain B cell activation due to receptor aggregation following binding of monovalent antigens presented on the antigen presenting cell.

Hierarchical graphs for rule-based modeling of biochemical systems

Background

In rule-based modeling, graphs are used to represent molecules: a colored vertex represents a component of a molecule, a vertex attribute represents the internal state of a component, and an edge represents a bond between components. Components of a molecule share the same color. Furthermore, graph-rewriting rules are used to represent molecular interactions. A rule that specifies addition (removal) of an edge represents a class of association (dissociation) reactions, and a rule that specifies a change of a vertex attribute represents a class of reactions that affect the internal state of a molecular component. A set of rules comprises an executable model that can be used to determine, through various means, the system-level dynamics of molecular interactions in a biochemical system.

Results

For purposes of model annotation, we propose the use of hierarchical graphs to represent structural relationships among components and subcomponents of molecules. We illustrate how hierarchical graphs can be used to naturally document the structural organization of the functional components and subcomponents of two proteins: the protein tyrosine kinase Lck and the T cell receptor (TCR) complex. We also show that computational methods developed for regular graphs can be applied to hierarchical graphs. In particular, we describe a generalization of Nauty, a graph isomorphism and canonical labeling algorithm. The generalized version of the Nauty procedure, which we call HNauty, can be used to assign canonical labels to hierarchical graphs or more generally to graphs with multiple edge types. The difference between the Nauty and HNauty procedures is minor, but for completeness, we provide an explanation of the entire HNauty algorithm.

Conclusions

Hierarchical graphs provide more intuitive formal representations of proteins and other structured molecules with multiple functional components than do the regular graphs of current languages for specifying rule-based models, such as the BioNetGen language (BNGL). Thus, the proposed use of hierarchical graphs should promote clarity and better understanding of rule-based models.

RuleMonkey: software for stochastic simulation of rule-based models

Background

The system-level dynamics of many molecular interactions, particularly protein-protein interactions, can be conveniently represented using reaction rules, which can be specified using model-specification languages, such as the BioNetGen language (BNGL). A set of rules implicitly defines a (bio)chemical reaction network. The reaction network implied by a set of rules is often very large, and as a result, generation of the network implied by rules tends to be computationally expensive. Moreover, the cost of many commonly used methods for simulating network dynamics is a function of network size. Together these factors have limited application of the rule-based modeling approach. Recently, several methods for simulating rule-based models have been developed that avoid the expensive step of network generation. The cost of these "network-free" simulation methods is independent of the number of reactions implied by rules. Software implementing such methods is now needed for the simulation and analysis of rule-based models of biochemical systems.

Results

Here, we present a software tool called RuleMonkey, which implements a network-free method for simulation of rule-based models that is similar to Gillespie's method. The method is suitable for rule-based models that can be encoded in BNGL, including models with rules that have global application conditions, such as rules for intramolecular association reactions. In addition, the method is rejection free, unlike other network-free methods that introduce null events, i.e., steps in the simulation procedure that do not change the state of the reaction system being simulated. We verify that RuleMonkey produces correct simulation results, and we compare its performance against DYNSTOC, another BNGL-compliant tool for network-free simulation of rule-based models. We also compare RuleMonkey against problem-specific codes implementing network-free simulation methods.

Conclusions

RuleMonkey enables the simulation of rule-based models for which the underlying reaction networks are large. It is typically faster than DYNSTOC for benchmark problems that we have examined. RuleMonkey is freely available as a stand-alone application http://public.tgen.org/rulemonkey. It is also available as a simulation engine within GetBonNie, a web-based environment for building, analyzing and sharing rule-based models.

Cannabinoid 1 receptor and interleukin-6 receptor together induce integration of protein kinase and transcription factor signaling to trigger neurite outgrowth

Activation of the G(o/i)-coupled cannabinoid 1 receptor (CB1R) has been shown to induce neurite outgrowth in Neuro2A cells through activation of Src kinase and STAT3 transcription factor. Signaling by the interleukin 6 receptor (IL-6R) also activates STAT3 through Jak kinase. We studied if signals from the two pathways could be integrated in a synergistic manner to trigger neurite outgrowth in Neuro2A cells. At low concentrations, when agonist at either receptor by itself has no effect, we found that CB1R and IL-6R stimulation together induced synergistic neurite outgrowth. Signal integration requires activation of transcription factors by Src, Jak, and mitogen-activated protein kinases. Mitogen-activated protein kinase can be activated by both receptors and shows enhanced early activation in the presence of both ligands. CREB and STAT3 transcription factors are required for synergy and show enhanced DNA-binding activity when both receptors are activated. STAT3 plays a critical role in integration of the signals downstream of the two receptors. When both pathways are activated, STAT3 phosphorylation is sustained for 6 h. This prolonged activation of STAT3 requires deactivation of SHP2 phosphatase. Reduction of SHP2 levels by RNA interference results in greater synergy in neurite outgrowth. Simultaneous knockdown of both SHP2 and STAT3 blocks the synergistic triggering of neurite outgrowth, indicating that STAT3 is downstream of SHP2. CB1R and IL-6R co-stimulation enhanced the differentiation of rat cortical neuron primary cultures. These results provide a mechanism where multiple protein kinases and transcription factors interact to integrate signals from G protein-coupled and cytokine receptor to evoke neurite outgrowth in Neuro2A cells.

A bipolar clamp mechanism for activation of Jak-family protein tyrosine kinases

Most cell surface receptors for growth factors and cytokines dimerize in order to mediate signal transduction. For many such receptors, the Janus kinase (Jak) family of non-receptor protein tyrosine kinases are recruited in pairs and juxtaposed by dimerized receptor complexes in order to activate one another by trans-phosphorylation. An alternative mechanism for Jak trans-phosphorylation has been proposed in which the phosphorylated kinase interacts with the Src homology 2 (SH2) domain of SH2-B, a unique adaptor protein with the capacity to homo-dimerize. Building on a rule-based kinetic modeling approach that considers the concerted nature and combinatorial complexity of modular protein domain interactions, we examine these mechanisms in detail, focusing on the growth hormone (GH) receptor/Jak2/SH2-Bβ system. The modeling results suggest that, whereas Jak2-(SH2-Bβ)₂-Jak2 heterotetramers are scarcely expected to affect Jak2 phosphorylation, SH2-Bβ and dimerized receptors synergistically promote Jak2 trans-activation in the context of intracellular signaling. Analysis of the results revealed a unique mechanism whereby SH2-B and receptor dimers constitute a bipolar ‘clamp’ that stabilizes the active configuration of two Jak2 molecules in the same macro-complex.

Janus kinases (Jaks) interact with and activate receptors on the cell surface that mediate changes in gene expression. How these interactions are promoted and regulated is of central interest in fields such as cellular endocrinology and immunology. Here, detailed computational models of Jak activation are offered at the level of protein modification states and interaction domains, wherein the specification of only a handful of binding/reaction rules can produce networks comprised of thousands of differential equations. Specifically, we investigated the role of an adaptor protein, SH2-B, revealing a novel mechanism whereby it cooperates with receptors to form a stable complex that juxtaposes two Jak molecules for efficient activation. We refer to this mode of molecular assembly as the bipolar clamp mechanism.

Complexity and modularity of intracellular networks: a systematic approach for modelling and simulation

Assembly of quantitative models of large complex networks brings about several challenges. One of them is the combinatorial complexity, where relatively few signalling molecules can combine to form thousands or millions of distinct chemical species. A receptor that has several separate phosphorylation sites can exist in hundreds of different states, many of which must be accounted for individually when simulating the time course of signalling. When assembly of protein complexes is being included, the number of distinct molecular species can easily increase by a few orders of magnitude. Validation, visualisation and understanding the network can become intractable. Another challenge appears when the modeller needs to recast or grow a model. Keeping track of changes and adding new elements present a significant difficulty. An approach to solve these challenges within the virtual cell (VCell) is described. Using (i) automatic extraction from pathway databases of model components (http://vcell.org/biopax) and (ii) rules of interactions that serve as reaction network generators (http://vcell.org/bionetgen), a way is provided for semi-automatic generation of quantitative mathematical models that also facilitates the reuse of model elements. In this approach, kinetic models of large, complex networks can be assembled from separately constructed modules, either directly or via rules. To implement this approach, the strength of several related technologies is combined: the BioPAX ontology, the BioNetGen rule-based description of molecular interactions and the VCell modelling and simulation framework.

Domain-oriented reduction of rule-based network models

The coupling of membrane-bound receptors to transcriptional regulators and other effector functions is mediated by multi-domain proteins that form complex assemblies. The modularity of protein interactions lends itself to a rule-based description, in which species and reactions are generated by rules that encode the necessary context for an interaction to occur, but also can produce a combinatorial explosion in the number of chemical species that make up the signalling network. The authors have shown previously that exact network reduction can be achieved using hierarchical control relationships between sites/domains on proteins to dissect multi-domain proteins into sets of non-interacting sites, allowing the replacement of each 'full' (progenitor) protein with a set of derived auxiliary (offspring) proteins. The description of a network in terms of auxiliary proteins that have fewer sites than progenitor proteins often greatly reduces network size. The authors describe here a method for automating domain-oriented model reduction and its implementation as a module in the BioNetGen modelling package. It takes as input a standard BioNetGen model and automatically performs the following steps: 1) detecting the hierarchical control relationships between sites; 2) building up the auxiliary proteins; 3) generating a raw reduced model and 4) cleaning up the raw model to provide the correct mass balance for each chemical species in the reduced network. The authors tested the performance of this module on models representing portions of growth factor receptor and immunoreceptor-mediated signalling networks and confirmed its ability to reduce the model size and simulation cost by at least one or two orders of magnitude. Limitations of the current algorithm include the inability to reduce models based on implicit site dependencies or heterodimerisation and loss of accuracy when dynamics are computed stochastically.

Rule-based modeling of biochemical systems with BioNetGen

Rule-based modeling involves the representation of molecules as structured objects and molecular interactions as rules for transforming the attributes of these objects. The approach is notable in that it allows one to systematically incorporate site-specific details about protein-protein interactions into a model for the dynamics of a signal-transduction system, but the method has other applications as well, such as following the fates of individual carbon atoms in metabolic reactions. The consequences of protein-protein interactions are difficult to specify and track with a conventional modeling approach because of the large number of protein phosphoforms and protein complexes that these interactions potentially generate. Here, we focus on how a rule-based model is specified in the BioNetGen language (BNGL) and how a model specification is analyzed using the BioNetGen software tool. We also discuss new developments in rule-based modeling that should enable the construction and analyses of comprehensive models for signal transduction pathways and similarly large-scale models for other biochemical systems.

Kinetic Monte Carlo method for rule-based modeling of biochemical networks

We present a kinetic Monte Carlo method for simulating chemical transformations specified by reaction rules, which can be viewed as generators of chemical reactions, or equivalently, definitions of reaction classes. A rule identifies the molecular components involved in a transformation, how these components change, conditions that affect whether a transformation occurs, and a rate law. The computational cost of the method, unlike conventional simulation approaches, is independent of the number of possible reactions, which need not be specified in advance or explicitly generated in a simulation. To demonstrate the method, we apply it to study the kinetics of multivalent ligand-receptor interactions. We expect the method will be useful for studying cellular signaling systems and other physical systems involving aggregation phenomena.

Computational models of tandem SRC homology 2 domain interactions and application to phosphoinositide 3-kinase

Intracellular signal transduction proteins typically utilize multiple interaction domains for proper targeting, and thus a broad diversity of distinct signaling complexes may be assembled. Considering the coordination of only two such domains, as in tandem Src homology 2 (SH2) domain constructs, gives rise to a kinetic scheme that is not adequately described by simple models used routinely to interpret in vitro binding measurements. To analyze the interactions between tandem SH2 domains and bisphosphorylated peptides, we formulated detailed kinetic models and applied them to the phosphoinositide 3-kinase p85 regulatory subunit/platelet-derived growth factor beta-receptor system. Data for this system from different in vitro assay platforms, including surface plasmon resonance, competition binding, and isothermal titration calorimetry, were reconciled to estimate the magnitude of the cooperativity characterizing the sequential binding of the high and low affinity SH2 domains (C-SH2 and N-SH2, respectively). Compared with values based on an effective volume approximation, the estimated cooperativity is 3 orders of magnitude lower, indicative of significant structural constraints. Homodimerization of full-length p85 was found to be an alternative mechanism for high avidity binding to phosphorylated platelet-derived growth factor receptors, which would render the N-SH2 domain dispensable for receptor binding.