Search tips
Search criteria

Results 1-25 (1016549)

Clipboard (0)

Related Articles

1.  Multiple Model-Informed Open-Loop Control of Uncertain Intracellular Signaling Dynamics 
PLoS Computational Biology  2014;10(4):e1003546.
Computational approaches to tune the activation of intracellular signal transduction pathways both predictably and selectively will enable researchers to explore and interrogate cell biology with unprecedented precision. Techniques to control complex nonlinear systems typically involve the application of control theory to a descriptive mathematical model. For cellular processes, however, measurement assays tend to be too time consuming for real-time feedback control and models offer rough approximations of the biological reality, thus limiting their utility when considered in isolation. We overcome these problems by combining nonlinear model predictive control with a novel adaptive weighting algorithm that blends predictions from multiple models to derive a compromise open-loop control sequence. The proposed strategy uses weight maps to inform the controller of the tendency for models to differ in their ability to accurately reproduce the system dynamics under different experimental perturbations (i.e. control inputs). These maps, which characterize the changing model likelihoods over the admissible control input space, are constructed using preexisting experimental data and used to produce a model-based open-loop control framework. In effect, the proposed method designs a sequence of control inputs that force the signaling dynamics along a predefined temporal response without measurement feedback while mitigating the effects of model uncertainty. We demonstrate this technique on the well-known Erk/MAPK signaling pathway in T cells. In silico assessment demonstrates that this approach successfully reduces target tracking error by 52% or better when compared with single model-based controllers and non-adaptive multiple model-based controllers. In vitro implementation of the proposed approach in Jurkat cells confirms a 63% reduction in tracking error when compared with the best of the single-model controllers. This study provides an experimentally-corroborated control methodology that utilizes the knowledge encoded within multiple mathematical models of intracellular signaling to design control inputs that effectively direct cell behavior in open-loop.
Author Summary
Most cell behavior arises as a response to external forces. Signals from the extracellular environment are passed to the cell's nucleus through a complex network of interacting proteins. Perturbing these pathways can change the strength or outcome of the signals, which could be used to treat or prevent a pathological response. While manipulating these networks can be achieved using a variety of methods, the ability to do so predictably over time would provide an unprecedented level of control over cell behavior and could lead to new therapeutic design and research tools in medicine and systems biology. Hence, we propose a practical computational framework to aid in the design of experimental perturbations to force cell signaling dynamics to follow a predefined response. Our approach represents a novel merger of model-based control and information theory to blend the predictions from multiple mathematical models into a meaningful compromise solution. We verify through simulation and experimentation that this solution produces excellent agreement between the cell readouts and several predefined trajectories, even in the presence of significant modeling uncertainty and without measurement feedback. By combining elements of information and control theory, our approach will help advance the best practices in model-based control applications for medicine.
PMCID: PMC3983080  PMID: 24722333
2.  A Canonical Model of Multistability and Scale-Invariance in Biological Systems 
PLoS Computational Biology  2012;8(8):e1002634.
Multistability and scale-invariant fluctuations occur in a wide variety of biological organisms from bacteria to humans as well as financial, chemical and complex physical systems. Multistability refers to noise driven switches between multiple weakly stable states. Scale-invariant fluctuations arise when there is an approximately constant ratio between the mean and standard deviation of a system's fluctuations. Both are an important property of human perception, movement, decision making and computation and they occur together in the human alpha rhythm, imparting it with complex dynamical behavior. Here, we elucidate their fundamental dynamical mechanisms in a canonical model of nonlinear bifurcations under stochastic fluctuations. We find that the co-occurrence of multistability and scale-invariant fluctuations mandates two important dynamical properties: Multistability arises in the presence of a subcritical Hopf bifurcation, which generates co-existing attractors, whilst the introduction of multiplicative (state-dependent) noise ensures that as the system jumps between these attractors, fluctuations remain in constant proportion to their mean and their temporal statistics become long-tailed. The simple algebraic construction of this model affords a systematic analysis of the contribution of stochastic and nonlinear processes to cortical rhythms, complementing a recently proposed biophysical model. Similar dynamics also occur in a kinetic model of gene regulation, suggesting universality across a broad class of biological phenomena.
Author Summary
Biological systems are able to adapt to rapidly and widely changing environments. Many biological organisms employ two distinct mechanisms that improve their survival in these circumstances: Firstly they exhibit rapid, qualitative changes in their internal dynamics; secondly they possess the ability to respond to change that is not absolute, but scales in proportion to the underlying intensity of the environment. In this paper, we study a simple class of noisy, dynamical systems that mathematically represent a very broad range of more complex models. We hence show how a combination of nonlinear instabilities and state-dependent noise in this model is able to unify these two apparently distinct biological phenomena. To illustrate its unifying potential, this simple model is applied to two very distinct biological processes – the spontaneous activity of the human cortex (i.e. when subjects are at rest), and genetic regulation in a bacteriophage. We also provide proof of principle that our model can be inverted from empirical data, allowing estimation of the parameters that express the nonlinear and stochastic influences at play in the underlying system.
PMCID: PMC3415415  PMID: 22912567
3.  A novel computational model of the circadian clock in Arabidopsis that incorporates PRR7 and PRR9 
We developed a mathematical model of the Arabidopsis circadian clock, including PRR7 and PRR9, which is able to predict several single, double and triple mutant phenotypes.Sensitivity Analysis was used to identify the properties and time sensing mechanisms of model structures.PRR7 and CCA1/LHY were identified as weak points of the mathematical model indicating where more experimental data is needed for further model development.Detailed dynamical studies showed that the timing of an evening light sensing element is essential for day length responsiveness
In recent years, molecular genetic techniques have revealed a complex network of components in the Arabidopsis circadian clock. Mathematical models allow for a detailed study of the dynamics and architecture of such complex gene networks leading to a better understanding of the genetic interactions. It is important to maintain a constant iteration with experimentation, to include novel components as they are discovered and use the updated model to design new experiments. This study develops a framework to introduce new components into the mathematical model of the Arabidopsis circadian clock accelerating the iterative model development process and gaining insight into the system's properties.
We used the interlocked feedback loop model published in Locke et al (2005) as the base model. In Arabidopsis, the first suggested regulatory loop involves the morning expressed transcription factors CIRCADIAN CLOCK-ASSOCIATED 1 (CCA1) and LATE ELONGATED HYPOCOTYL (LHY), and the evening expressed pseudo-response regulator TIMING OF CAB EXPRESSION (TOC1). The hypothetical component X had been introduced to realize a longer delay between gene expression of CCA1/LHY and TOC1. The introduction of Y was motivated by the need for a mechanism to reproduce the dampening short period rhythms of the cca1/lhy double mutant and to include an additional light input at the end of the day.
In this study, the new components pseudo-response regulators PRR7 and PRR9 were added in negative feedback loops based on the biological hypothesis that they are activated by LHY and in turn repress LHY transcription (Farré et al, 2005; Figure 1). We present three iterations steps of model development (Figure 1A–C).
A wide range of tools was used to establish and analyze new model structures. One of the challenges facing mathematical modeling of biological processes is parameter identification; they are notoriously difficult to determine experimentally. We established an optimization procedure based on an evolutionary strategy with a cost function mainly derived from wild-type characteristics. This ensured that the model was not restricted by a specific set of parameters and enabled us to use a large set of biological mutant information to assess the predictive capability of the model structure. Models were evaluated by means of an extended phenotype catalogue, allowing for an easy and fair comparison of the structures. We also carried out detailed simulation analysis of component interactions to identify weak points in the structure and suggest further modifications. Finally, we applied sensitivity analysis in a novel manner, using it to direct the model development. Sensitivity analysis provides quantitative measures of robustness; the two measures in this study were the traces of component concentrations over time (classical state sensitivities) and phase behavior (measured by the phase response curve). Three major results emerged from the model development process.
First, the iteration process helped us to learn about general characteristics of the system. We observed that the timing of Y expression is critical for evening light entrainment, which enables the system to respond to changes in day length. This is important for our understanding of the mechanism of light input to the clock and will add in the identification of biological candidates for this function. In addition, our results suggest that a detailed description of the mechanisms of genetic interactions is important for the systems behavior. We observed that the introduction of an experimentally based precise light regulation mechanism on PRR9 expression had a significant effect on the systems behavior.
Second, the final model structure (Figure 1C) was capable of predicting a wide range of mutant phenotypes, such as a reduction of TOC1 expression by RNAi (toc1RNAi), mutations in PRR7 and PRR9 and the novel mutant combinations prr9toc1RNAi and prr7prr9toc1RNAi. However, it was unable to predict the mutations in CCA1 and LHY.
Finally, sensitivity analysis identified the weak points of the system. The developed model structure was heavily based on the TOC1/Y feedback loop. This could explain the model's failure to predict the cca1lhy double mutant phenotype. More detailed information on the regulation of CCA1 and LHY expression will be important to achieve the right balance between the different regulatory loops in the mathematical model. This is in accordance with genetic studies that have identified several genes involved in the regulation of LHY and CCA1 expression. The identification of their mechanism of action will be necessary for the next model development.
In plants, as in animals, the core mechanism to retain rhythmic gene expression relies on the interaction of multiple feedback loops. In recent years, molecular genetic techniques have revealed a complex network of clock components in Arabidopsis. To gain insight into the dynamics of these interactions, new components need to be integrated into the mathematical model of the plant clock. Our approach accelerates the iterative process of model identification, to incorporate new components, and to systematically test different proposed structural hypotheses. Recent studies indicate that the pseudo-response regulators PRR7 and PRR9 play a key role in the core clock of Arabidopsis. We incorporate PRR7 and PRR9 into an existing model involving the transcription factors TIMING OF CAB (TOC1), LATE ELONGATED HYPOCOTYL (LHY) and CIRCADIAN CLOCK ASSOCIATED (CCA1). We propose candidate models based on experimental hypotheses and identify the computational models with the application of an optimization routine. Validation is accomplished through systematic analysis of various mutant phenotypes. We introduce and apply sensitivity analysis as a novel tool for analyzing and distinguishing the characteristics of proposed architectures, which also allows for further validation of the hypothesized structures.
PMCID: PMC1682023  PMID: 17102803
Arabidopsis; circadian rhythms; mathematical modeling; parameter optimization; sensitivity analysis
4.  Dynamic optimization of distributed biological systems using robust and efficient numerical techniques 
BMC Systems Biology  2012;6:79.
Systems biology allows the analysis of biological systems behavior under different conditions through in silico experimentation. The possibility of perturbing biological systems in different manners calls for the design of perturbations to achieve particular goals. Examples would include, the design of a chemical stimulation to maximize the amplitude of a given cellular signal or to achieve a desired pattern in pattern formation systems, etc. Such design problems can be mathematically formulated as dynamic optimization problems which are particularly challenging when the system is described by partial differential equations.
This work addresses the numerical solution of such dynamic optimization problems for spatially distributed biological systems. The usual nonlinear and large scale nature of the mathematical models related to this class of systems and the presence of constraints on the optimization problems, impose a number of difficulties, such as the presence of suboptimal solutions, which call for robust and efficient numerical techniques.
Here, the use of a control vector parameterization approach combined with efficient and robust hybrid global optimization methods and a reduced order model methodology is proposed. The capabilities of this strategy are illustrated considering the solution of a two challenging problems: bacterial chemotaxis and the FitzHugh-Nagumo model.
In the process of chemotaxis the objective was to efficiently compute the time-varying optimal concentration of chemotractant in one of the spatial boundaries in order to achieve predefined cell distribution profiles. Results are in agreement with those previously published in the literature. The FitzHugh-Nagumo problem is also efficiently solved and it illustrates very well how dynamic optimization may be used to force a system to evolve from an undesired to a desired pattern with a reduced number of actuators. The presented methodology can be used for the efficient dynamic optimization of generic distributed biological systems.
PMCID: PMC3532319  PMID: 22748139
Dynamic optimization; Distributed biological systems; Reduced order models; Global optimization methods; Hybrid optimization methods; Pattern formation and control
5.  Nonlinear systems in medicine. 
Many achievements in medicine have come from applying linear theory to problems. Most current methods of data analysis use linear models, which are based on proportionality between two variables and/or relationships described by linear differential equations. However, nonlinear behavior commonly occurs within human systems due to their complex dynamic nature; this cannot be described adequately by linear models. Nonlinear thinking has grown among physiologists and physicians over the past century, and non-linear system theories are beginning to be applied to assist in interpreting, explaining, and predicting biological phenomena. Chaos theory describes elements manifesting behavior that is extremely sensitive to initial conditions, does not repeat itself and yet is deterministic. Complexity theory goes one step beyond chaos and is attempting to explain complex behavior that emerges within dynamic nonlinear systems. Nonlinear modeling still has not been able to explain all of the complexity present in human systems, and further models still need to be refined and developed. However, nonlinear modeling is helping to explain some system behaviors that linear systems cannot and thus will augment our understanding of the nature of complex dynamic systems within the human body in health and in disease states.
PMCID: PMC2588816  PMID: 14580107
6.  Quantitative analysis of regulatory flexibility under changing environmental conditions 
Day length changes with the seasons in temperate latitudes, affecting the many biological rhythms that entrain to the day/night cycle: we measure these effects on the expression of Arabidopsis clock genes, using RNA and reporter gene readouts, with a new method of phase analysis.Dusk sensitivity is proposed as a simple, natural and general mathematical measure to analyse and manipulate the changing phase of a clock output relative to the change in the day/night cycle.Dusk sensitivity shows how increasing the numbers of feedback loops in the Arabidopsis clock models allows more flexible regulation, consistent with a previously-proposed, general operating principle of biological networks.The Arabidopsis clock genes show flexibility of regulation that is characteristic of a three-loop clock model, validating aspects of the model and the operating principle, but some clock output genes show greater flexibility arising from direct light regulation.
The analysis of dynamic, non-linear regulation with the aid of mechanistic models is central to Systems Biology. This study compares the predictions of mechanistic, mathematical models of the circadian clock with molecular time-series data on rhythmic gene expression in the higher plant Arabidopsis thaliana. Analysis of the models helps us to understand (explain and predict) how the clock gene circuit balances regulation by external and endogenous factors to achieve particular behaviours. Such multi-factorial regulation is ubiquitous in, and characteristic of, living systems.
The Earth's rotation causes predictable changes in the environment, notably in the availability of sunlight for photosynthesis. Many biological processes are driven by the environmental input via sensory pathways, for example, from photoreceptors. Circadian clocks provide an alternative strategy. These endogenous, 24-h rhythms can drive biological processes that anticipate the regular environmental changes, rather than merely responding. Many rhythmic processes have both light and clock control. Indeed, the clock components themselves must balance internal timing with external inputs, because circadian clocks are reset daily through light regulation of one or more clock components. This process of entrainment is complicated by the change in day length. When the times of dawn and dusk move apart in summer, and closer together in winter, does the clock track dawn, track dusk or interpolate between them?
In plants, the clock controls leaf and petal movements, the opening and closing of stomatal pores, the discharge of floral fragrances, and many metabolic activities, especially those associated with photosynthesis. Centuries of physiological studies have shown that these rhythms can behave differently. Flowering in Ipomoea nil (Pharbitis nil, Japanese morning glory) is controlled by a rhythm that tracks the time of dusk, to give a classic example. We showed that two other rhythms associated with vegetative growth track dawn in this species (Figure 5A), so the clock system allows flexible regulation.
The relatively small number of components involved in the circadian clockwork makes it an ideal candidate for mathematical modelling. Molecular genetic studies in a variety of model eukaryotes have shown that the circadian rhythm is generated by a network of 6–20 genes. These genes form feedback loops generating a rhythm in mRNA production. A single negative feedback loop in which a gene encodes a protein that, after several hours, turns off transcription is capable of generating a circadian rhythm, in principle. A single light input can entrain the clock to ‘local time', synchronised with a light–dark cycle. However, real circadian clocks have proven to be more complicated than this, with multiple light inputs and interlocked feedback loops.
We have previously argued from mathematical analysis that multi-loop networks increase the flexibility of regulation (Rand et al, 2004) and have shown that appropriately deployed flexibility can confer functional robustness (Akman et al, 2010). Here we test whether that flexibility can be demonstrated in vivo, in the model plant, A. thaliana. The Arabidopsis clock mechanism comprises a feedback loop in which two partially redundant, myb transcription factors, LATE ELONGATED HYPOCOTYL (LHY) and CIRCADIAN CLOCK ASSOCIATED 1 (CCA1), repress the expression of their activator, TIMING OF CAB EXPRESSION 1 (TOC1). We previously modelled this single-loop circuit and showed that it was not capable of recreating important data (Locke et al, 2005a). An extended, two-loop model was developed to match observed behaviours, incorporating a hypothetical gene Y, for which the best identified candidate was the GIGANTEA gene (GI) (Locke et al, 2005b). Two further models incorporated the TOC1 homologues PSEUDO-RESPONSE REGULATOR (PRR) 9 and PRR7 (Locke et al, 2006; Zeilinger et al, 2006). In these circuits, a morning oscillator (LHY/CCA1–PRR9/7) is coupled to an evening oscillator (Y/GI–TOC1) via the original LHY/CCA1–TOC1 loop.
These clock models, like those for all other organisms, were developed using data from simple conditions of constant light, darkness or 12-h light–12-h dark cycles. We therefore tested how the clock genes in Arabidopsis responded to light–dark cycles with different photoperiods, from 3 h light to 18 h light per 24-h cycle (Edinburgh, 56° North latitude, has 17.5 h light in midsummer). The time-series assays of mRNA and in vivo reporter gene images showed a range of peak times for different genes, depending on the photoperiod (Figure 5C). A new data analysis method, mFourfit, was introduced to measure the peak times, in the Biological Rhythms Analysis Software Suite (BRASS v3.0). None of the genes showed the dusk-tracking behaviour characteristic of the Ipomoea flowering rhythm. The one-, two- and three-loop models were analysed to understand the observed patterns. A new mathematical measure, dusk sensitivity, was introduced to measure the change in timing of a model component versus a change in the time of dusk. The one- and two-loop models tracked dawn and dusk, respectively, under all conditions. Only the three-loop model (Figure 5B) had the flexibility required to match the photoperiod-dependent changes that we found in vivo, and in particular the unexpected, V-shaped pattern in the peak time of TOC1 expression. This pattern of regulation depends on the structure and light inputs to the model's evening oscillator, so the in vivo data supported this aspect of the model. LHY and CCA1 gene expression under short photoperiods showed greater dusk sensitivity, in the interval 2–6 h before dawn, than the three-loop model predicted, so these data will help to constrain future models.
The approach described here could act as a template for experimental biologists seeking to understand biological regulation using dynamic, experimental perturbations and time-series data. Simulation of mathematical models (despite known imperfections) can provide contrasting hypotheses that guide understanding. The system's detailed behaviour is complex, so a natural and general measure such as dusk sensitivity is helpful to focus on one property of the system. We used the measure to compare models, and to predict how this property could be manipulated. To enable additional analysis of this system, we provide the time-series data and experimental metadata online.
The circadian clock controls 24-h rhythms in many biological processes, allowing appropriate timing of biological rhythms relative to dawn and dusk. Known clock circuits include multiple, interlocked feedback loops. Theory suggested that multiple loops contribute the flexibility for molecular rhythms to track multiple phases of the external cycle. Clear dawn- and dusk-tracking rhythms illustrate the flexibility of timing in Ipomoea nil. Molecular clock components in Arabidopsis thaliana showed complex, photoperiod-dependent regulation, which was analysed by comparison with three contrasting models. A simple, quantitative measure, Dusk Sensitivity, was introduced to compare the behaviour of clock models with varying loop complexity. Evening-expressed clock genes showed photoperiod-dependent dusk sensitivity, as predicted by the three-loop model, whereas the one- and two-loop models tracked dawn and dusk, respectively. Output genes for starch degradation achieved dusk-tracking expression through light regulation, rather than a dusk-tracking rhythm. Model analysis predicted which biochemical processes could be manipulated to extend dusk tracking. Our results reveal how an operating principle of biological regulators applies specifically to the plant circadian clock.
PMCID: PMC3010117  PMID: 21045818
Arabidopsis thaliana; biological clocks; dynamical systems; gene regulatory networks; mathematical models; photoperiodism
7.  The Origins of Time-Delay in Template Biopolymerization Processes 
PLoS Computational Biology  2010;6(4):e1000726.
Time-delays are common in many physical and biological systems and they give rise to complex dynamic phenomena. The elementary processes involved in template biopolymerization, such as mRNA and protein synthesis, introduce significant time delays. However, there is not currently a systematic mapping between the individual mechanistic parameters and the time delays in these networks. We present here the development of mathematical, time-delay models for protein translation, based on PDE models, which in turn are derived through systematic approximations of first-principles mechanistic models. Theoretical analysis suggests that the key features that determine the time-delays and the agreement between the time-delay and the mechanistic models are ribosome density and distribution, i.e., the number of ribosomes on the mRNA chain relative to their maximum and their distribution along the mRNA chain. Based on analytical considerations and on computational studies, we show that the steady-state and dynamic responses of the time-delay models are in excellent agreement with the detailed mechanistic models, under physiological conditions that correspond to uniform ribosome distribution and for ribosome density up to 70%. The methodology presented here can be used for the development of reduced time-delay models of mRNA synthesis and large genetic networks. The good agreement between the time-delay and the mechanistic models will allow us to use the reduced model and advanced computational methods from nonlinear dynamics in order to perform studies that are not practical using the large-scale mechanistic models.
Author Summary
Genetic networks display exceedingly complex and rich behavior which is modulated by multiple mechanisms, including many diverse types of interactions between DNA, mRNA and protein molecules. Mathematical models of gene networks must necessarily consider the essential mechanistic details of the processes involved in order to make reliable predictions. However, even though the description of the process becomes more accurate as more mechanistic details are incorporated into the mathematical model, the added mathematical complexity will make it difficult to parameterize and extract information from such models given the limited amount of experimental data. Protein synthesis is precisely one of the phases in the network machinery where certain mechanistic details are important and should thus be taken into account. Here, we develop a methodology to reduce a mathematical model for protein synthesis by performing approximations on a mechanistic model, retaining the essential details of the process. Our methodology opens up the possibility of utilizing powerful mathematical tools, such as bifurcation analysis, for understanding the complex dynamics displayed by genetic networks and design strategies for metabolic engineering and synthetic biology.
PMCID: PMC2848540  PMID: 20369012
8.  Modeling Integrated Cellular Machinery Using Hybrid Petri-Boolean Networks 
PLoS Computational Biology  2013;9(11):e1003306.
The behavior and phenotypic changes of cells are governed by a cellular circuitry that represents a set of biochemical reactions. Based on biological functions, this circuitry is divided into three types of networks, each encoding for a major biological process: signal transduction, transcription regulation, and metabolism. This division has generally enabled taming computational complexity dealing with the entire system, allowed for using modeling techniques that are specific to each of the components, and achieved separation of the different time scales at which reactions in each of the three networks occur. Nonetheless, with this division comes loss of information and power needed to elucidate certain cellular phenomena. Within the cell, these three types of networks work in tandem, and each produces signals and/or substances that are used by the others to process information and operate normally. Therefore, computational techniques for modeling integrated cellular machinery are needed. In this work, we propose an integrated hybrid model (IHM) that combines Petri nets and Boolean networks to model integrated cellular networks. Coupled with a stochastic simulation mechanism, the model simulates the dynamics of the integrated network, and can be perturbed to generate testable hypotheses. Our model is qualitative and is mostly built upon knowledge from the literature and requires fine-tuning of very few parameters. We validated our model on two systems: the transcriptional regulation of glucose metabolism in human cells, and cellular osmoregulation in S. cerevisiae. The model produced results that are in very good agreement with experimental data, and produces valid hypotheses. The abstract nature of our model and the ease of its construction makes it a very good candidate for modeling integrated networks from qualitative data. The results it produces can guide the practitioner to zoom into components and interconnections and investigate them using such more detailed mathematical models.
Author Summary
Within the cell of an organism, three networks—signaling, transcriptional, and metabolic—are always at work to determine the response of the cell to signals from its environment, and consequently, its fate. Evidence from experimental studies is painting a picture of complex crosstalk among these networks. Thus, while a wide array of computational techniques exist for analyzing each of these network types, there is clear need for new modeling techniques that allow for simultaneously analyzing integrated networks, which combine elements from all three networks. Here, we provide a step towards achieving this task by combining two population modeling techniques—Petri nets and Boolean networks—to produce an integrated hybrid model. We demonstrate the accuracy and utility of this model on two biological systems: transcriptional regulation of glucose metabolism in human cells, and cellular osmoregulation in yeast.
PMCID: PMC3820535  PMID: 24244124
9.  Lineage grammars: describing, simulating and analyzing population dynamics 
BMC Bioinformatics  2014;15:249.
Precise description of the dynamics of biological processes would enable the mathematical analysis and computational simulation of complex biological phenomena. Languages such as Chemical Reaction Networks and Process Algebras cater for the detailed description of interactions among individuals and for the simulation and analysis of ensuing behaviors of populations. However, often knowledge of such interactions is lacking or not available. Yet complete oblivion to the environment would make the description of any biological process vacuous. Here we present a language for describing population dynamics that abstracts away detailed interaction among individuals, yet captures in broad terms the effect of the changing environment, based on environment-dependent Stochastic Tree Grammars (eSTG). It is comprised of a set of stochastic tree grammar transition rules, which are context-free and as such abstract away specific interactions among individuals. Transition rule probabilities and rates, however, can depend on global parameters such as population size, generation count, and elapsed time.
We show that eSTGs conveniently describe population dynamics at multiple levels including cellular dynamics, tissue development and niches of organisms. Notably, we show the utilization of eSTG for cases in which the dynamics is regulated by environmental factors, which affect the fate and rate of decisions of the different species. eSTGs are lineage grammars, in the sense that execution of an eSTG program generates the corresponding lineage trees, which can be used to analyze the evolutionary and developmental history of the biological system under investigation. These lineage trees contain a representation of the entire events history of the system, including the dynamics that led to the existing as well as to the extinct individuals.
We conclude that our suggested formalism can be used to easily specify, simulate and analyze complex biological systems, and supports modular description of local biological dynamics that can be later used as “black boxes” in a larger scope, thus enabling a gradual and hierarchical definition and simulation of complex biological systems. The simple, yet robust formalism enables to target a broad class of stochastic dynamic behaviors, especially those that can be modeled using global environmental feedback regulation rather than direct interaction between individuals.
PMCID: PMC4223406  PMID: 25047682
10.  Quantifying Transient 3D Dynamical Phenomena of Single mRNA Particles in Live Yeast Cell Measurements 
The journal of physical chemistry. B  2013;117(49):10.1021/jp4064214.
Single-particle tracking (SPT) has been extensively used to obtain information about diffusion and directed motion in a wide range of biological applications. Recently, new methods have appeared for obtaining precise (10s of nm) spatial information in three dimensions (3D) with high temporal resolution (measurements obtained every 4ms), which promise to more accurately sense the true dynamical behavior in the natural 3D cellular environment. Despite the quantitative 3D tracking information, the range of mathematical methods for extracting information about the underlying system has been limited mostly to mean-squared displacement analysis and other techniques not accounting for complex 3D kinetic interactions. There is a great need for new analysis tools aiming to more fully extract the biological information content from in vivo SPT measurements. High-resolution SPT experimental data has enormous potential to objectively scrutinize various proposed mechanistic schemes arising from theoretical biophysics and cell biology. At the same time, methods for rigorously checking the statistical consistency of both model assumptions and estimated parameters against observed experimental data (i.e. goodness-of-fit tests) have not received great attention. We demonstrate methods enabling (1) estimation of the parameters of 3D stochastic differential equation (SDE) models of the underlying dynamics given only one trajectory; and (2) construction of hypothesis tests checking the consistency of the fitted model with the observed trajectory so that extracted parameters are not over-interpreted (the tools are applicable to linear or nonlinear SDEs calibrated from non-stationary time series data). The approach is demonstrated on high-resolution 3D trajectories of single ARG3 mRNA particles in yeast cells in order to show the power of the methods in detecting signatures of transient directed transport. The methods presented are generally relevant to a wide variety of 2D and 3D SPT tracking applications.
PMCID: PMC3865222  PMID: 24015725
single particle tracking; stochastic differential equations; diffusion; goodness-of-fit testing; 3D microscopy
11.  Hierarchical Cluster-based Partial Least Squares Regression (HC-PLSR) is an efficient tool for metamodelling of nonlinear dynamic models 
BMC Systems Biology  2011;5:90.
Deterministic dynamic models of complex biological systems contain a large number of parameters and state variables, related through nonlinear differential equations with various types of feedback. A metamodel of such a dynamic model is a statistical approximation model that maps variation in parameters and initial conditions (inputs) to variation in features of the trajectories of the state variables (outputs) throughout the entire biologically relevant input space. A sufficiently accurate mapping can be exploited both instrumentally and epistemically. Multivariate regression methodology is a commonly used approach for emulating dynamic models. However, when the input-output relations are highly nonlinear or non-monotone, a standard linear regression approach is prone to give suboptimal results. We therefore hypothesised that a more accurate mapping can be obtained by locally linear or locally polynomial regression. We present here a new method for local regression modelling, Hierarchical Cluster-based PLS regression (HC-PLSR), where fuzzy C-means clustering is used to separate the data set into parts according to the structure of the response surface. We compare the metamodelling performance of HC-PLSR with polynomial partial least squares regression (PLSR) and ordinary least squares (OLS) regression on various systems: six different gene regulatory network models with various types of feedback, a deterministic mathematical model of the mammalian circadian clock and a model of the mouse ventricular myocyte function.
Our results indicate that multivariate regression is well suited for emulating dynamic models in systems biology. The hierarchical approach turned out to be superior to both polynomial PLSR and OLS regression in all three test cases. The advantage, in terms of explained variance and prediction accuracy, was largest in systems with highly nonlinear functional relationships and in systems with positive feedback loops.
HC-PLSR is a promising approach for metamodelling in systems biology, especially for highly nonlinear or non-monotone parameter to phenotype maps. The algorithm can be flexibly adjusted to suit the complexity of the dynamic model behaviour, inviting automation in the metamodelling of complex systems.
PMCID: PMC3127793  PMID: 21627852
12.  Identification of New IκBα Complexes by an Iterative Experimental and Mathematical Modeling Approach 
PLoS Computational Biology  2014;10(3):e1003528.
The transcription factor nuclear factor kappa-B (NFκB) is a key regulator of pro-inflammatory and pro-proliferative processes. Accordingly, uncontrolled NFκB activity may contribute to the development of severe diseases when the regulatory system is impaired. Since NFκB can be triggered by a huge variety of inflammatory, pro-and anti-apoptotic stimuli, its activation underlies a complex and tightly regulated signaling network that also includes multi-layered negative feedback mechanisms. Detailed understanding of this complex signaling network is mandatory to identify sensitive parameters that may serve as targets for therapeutic interventions. While many details about canonical and non-canonical NFκB activation have been investigated, less is known about cellular IκBα pools that may tune the cellular NFκB levels. IκBα has so far exclusively been described to exist in two different forms within the cell: stably bound to NFκB or, very transiently, as unbound protein. We created a detailed mathematical model to quantitatively capture and analyze the time-resolved network behavior. By iterative refinement with numerous biological experiments, we yielded a highly identifiable model with superior predictive power which led to the hypothesis of an NFκB-lacking IκBα complex that contains stabilizing IKK subunits. We provide evidence that other but canonical pathways exist that may affect the cellular IκBα status. This additional IκBα:IKKγ complex revealed may serve as storage for the inhibitor to antagonize undesired NFκB activation under physiological and pathophysiological conditions.
Author Summary
In unstimulated cells, the transcription factor NFκB resides in the cytosol bound to its inhibitor IκBα. Canonical activation of NFκB by numerous stimuli leads to proteasomal depletion of IκBα, thereby liberating NFκB to translocate into the nucleus to induce transcription of genes leading to proliferation, angiogenesis, metastasis, or chronic inflammation. Consequently, only transient activity needs to be warranted by immediate NFκB-dependent induction of negative regulatory mechanisms, including up-regulation of its inhibitor IκBα. Resynthesized IκBα consequently terminates NFκB activity by binding to its nuclear localization sequence. However, under physiological or pathophysiological conditions, random NFκB activation may occur, which needs to be avoided in order to guarantee proper cellular function. Using detailed dynamical modeling, we have now identified an additional IκBα containing complex to exist in un-stimulated cells which lacks NFκB but includes IKKγ (IκBα:IKKγ complex). This additional IκBα is not depleted from cells in the canonical fashion and may therefore serve as a cellular backup to avoid random NFκB activation.
PMCID: PMC3967930  PMID: 24675998
13.  Characterization of an inducible promoter in different DNA copy number conditions 
BMC Bioinformatics  2012;13(Suppl 4):S11.
The bottom-up programming of living organisms to implement novel user-defined biological capabilities is one of the main goals of synthetic biology. Currently, a predominant problem connected with the construction of even simple synthetic biological systems is the unpredictability of the genetic circuitry when assembled and incorporated in living cells. Copy number, transcriptional/translational demand and toxicity of the DNA-encoded functions are some of the major factors which may lead to cell overburdening and thus to nonlinear effects on system output. It is important to disclose the linearity working boundaries of engineered biological systems when dealing with such phenomena.
The output of an N-3-oxohexanoyl-L-homoserine lactone (HSL)-inducible RFP-expressing device was studied in Escherichia coli in different copy number contexts, ranging from 1 copy per cell (integrated in the genome) to hundreds (via multicopy plasmids). The system is composed by a luxR constitutive expression cassette and a RFP gene regulated by the luxI promoter, which is activated by the HSL-LuxR complex. System output, in terms of promoter activity as a function of HSL concentration, was assessed relative to the one of a reference promoter in identical conditions by using the Relative Promoter Units (RPU) approach. Nonlinear effects were observed in the maximum activity, which is identical in single and low copy conditions, while it decreases for higher copy number conditions. In order to properly compare the luxI promoter strength among all the conditions, a mathematical modeling approach was used to relate the promoter activity to the estimated HSL-LuxR complex concentration, which is the actual activator of transcription. During model fitting, a correlation between the copy number and the dissociation constant of HSL-LuxR complex and luxI promoter was observed.
Even in a simple inducible system, nonlinear effects are observed and non-trivial data processing is necessary to fully characterize its operation. The in-depth analysis of model systems like this can contribute to the advances in the synthetic biology field, since increasing the knowledge about linearity and working boundaries of biological phenomena could lead to a more rational design of artificial systems, also through mathematical models, which, for example, have been used here to study hard-to-predict interactions.
PMCID: PMC3314568  PMID: 22536957
14.  Analytical Derivation of Moment Equations in Stochastic Chemical Kinetics 
Chemical engineering science  2011;66(3):268-277.
The master probability equation captures the dynamic behavior of a variety of stochastic phenomena that can be modeled as Markov processes. Analytical solutions to the master equation are hard to come by though because they require the enumeration of all possible states and the determination of the transition probabilities between any two states. These two tasks quickly become intractable for all but the simplest of systems. Instead of determining how the probability distribution changes in time, we can express the master probability distribution as a function of its moments, and, we can then write transient equations for the probability distribution moments. In 1949, Moyal defined the derivative, or jump, moments of the master probability distribution. These are measures of the rate of change in the probability distribution moment values, i.e. what the impact is of any given transition between states on the moment values. In this paper we present a general scheme for deriving analytical moment equations for any N-dimensional Markov process as a function of the jump moments. Importantly, we propose a scheme to derive analytical expressions for the jump moments for any N-dimensional Markov process. To better illustrate the concepts, we focus on stochastic chemical kinetics models for which we derive analytical relations for jump moments of arbitrary order. Chemical kinetics models are widely used to capture the dynamic behavior of biological systems. The elements in the jump moment expressions are a function of the stoichiometric matrix and the reaction propensities, i.e the probabilistic reaction rates. We use two toy examples, a linear and a non-linear set of reactions, to demonstrate the applicability and limitations of the scheme. Finally, we provide an estimate on the minimum number of moments necessary to obtain statistical significant data that would uniquely determine the dynamics of the underlying stochastic chemical kinetic system. The first two moments only provide limited information, especially when complex, non-linear dynamics are involved.
PMCID: PMC3176737  PMID: 21949443
Multiscale models; Master equation; Model reduction; Statistical thermodynamics; Mathematical modelling; Computational chemistry
15.  Robust synchronization control scheme of a population of nonlinear stochastic synthetic genetic oscillators under intrinsic and extrinsic molecular noise via quorum sensing 
BMC Systems Biology  2012;6:136.
Collective rhythms of gene regulatory networks have been a subject of considerable interest for biologists and theoreticians, in particular the synchronization of dynamic cells mediated by intercellular communication. Synchronization of a population of synthetic genetic oscillators is an important design in practical applications, because such a population distributed over different host cells needs to exploit molecular phenomena simultaneously in order to emerge a biological phenomenon. However, this synchronization may be corrupted by intrinsic kinetic parameter fluctuations and extrinsic environmental molecular noise. Therefore, robust synchronization is an important design topic in nonlinear stochastic coupled synthetic genetic oscillators with intrinsic kinetic parameter fluctuations and extrinsic molecular noise.
Initially, the condition for robust synchronization of synthetic genetic oscillators was derived based on Hamilton Jacobi inequality (HJI). We found that if the synchronization robustness can confer enough intrinsic robustness to tolerate intrinsic parameter fluctuation and extrinsic robustness to filter the environmental noise, then robust synchronization of coupled synthetic genetic oscillators is guaranteed. If the synchronization robustness of a population of nonlinear stochastic coupled synthetic genetic oscillators distributed over different host cells could not be maintained, then robust synchronization could be enhanced by external control input through quorum sensing molecules. In order to simplify the analysis and design of robust synchronization of nonlinear stochastic synthetic genetic oscillators, the fuzzy interpolation method was employed to interpolate several local linear stochastic coupled systems to approximate the nonlinear stochastic coupled system so that the HJI-based synchronization design problem could be replaced by a simple linear matrix inequality (LMI)-based design problem, which could be solved with the help of LMI toolbox in MATLAB easily.
If the synchronization robustness criterion, i.e. the synchronization robustness ≥ intrinsic robustness + extrinsic robustness, then the stochastic coupled synthetic oscillators can be robustly synchronized in spite of intrinsic parameter fluctuation and extrinsic noise. If the synchronization robustness criterion is violated, external control scheme by adding inducer can be designed to improve synchronization robustness of coupled synthetic genetic oscillators. The investigated robust synchronization criteria and proposed external control method are useful for a population of coupled synthetic networks with emergent synchronization behavior, especially for multi-cellular, engineered networks.
PMCID: PMC3554485  PMID: 23101662
16.  Computational Cellular Dynamics Based on the Chemical Master Equation: A Challenge for Understanding Complexity 
Modern molecular biology has always been a great source of inspiration for computational science. Half a century ago, the challenge from understanding macromolecular dynamics has led the way for computations to be part of the tool set to study molecular biology. Twenty-five years ago, the demand from genome science has inspired an entire generation of computer scientists with an interest in discrete mathematics to join the field that is now called bioinformatics. In this paper, we shall lay out a new mathematical theory for dynamics of biochemical reaction systems in a small volume (i.e., mesoscopic) in terms of a stochastic, discrete-state continuous-time formulation, called the chemical master equation (CME). Similar to the wavefunction in quantum mechanics, the dynamically changing probability landscape associated with the state space provides a fundamental characterization of the biochemical reaction system. The stochastic trajectories of the dynamics are best known through the simulations using the Gillespie algorithm. In contrast to the Metropolis algorithm, this Monte Carlo sampling technique does not follow a process with detailed balance. We shall show several examples how CMEs are used to model cellular biochemical systems. We shall also illustrate the computational challenges involved: multiscale phenomena, the interplay between stochasticity and nonlinearity, and how macroscopic determinism arises from mesoscopic dynamics. We point out recent advances in computing solutions to the CME, including exact solution of the steady state landscape and stochastic differential equations that offer alternatives to the Gilespie algorithm. We argue that the CME is an ideal system from which one can learn to understand “complex behavior” and complexity theory, and from which important biological insight can be gained.
PMCID: PMC4079062  PMID: 24999297
biochemical networks; cellular signaling; epigenetics; master equation; nonlinear reactions; stochastic modeling
17.  Synthetic in vitro transcriptional oscillators 
A fundamental goal of synthetic biology is to understand design principles through engineering biochemical systems.Three in vitro synthetic transcriptional oscillators were constructed and analyzed: a two-node-negative feedback oscillator, an amplified negative-feedback oscillator, and a three-node ring oscillator.The in vitro oscillators are governed by similar design principles as previous theoretical studies and synthetic oscillators in vivo.Because of unintended reactions that arise even without the complexity of living cells, several challenges remain for predictive and robust oscillator performance.
Fundamental goals for synthetic biology are to understand the principles of biological circuitry from an engineering perspective and to establish engineering methods for creating biochemical circuitry to control molecular processes—both in vitro and in vivo (Benner and Sismour, 2005; Adrianantoandro et al, 2006). Here, we make use of a previously proposed class of in vitro biochemical systems, transcriptional circuits, that can be modularly wired into arbitrarily complex networks by changing the regulatory and coding sequence domains of DNA templates (Kim et al, 2006; Subsoontorn et al 2011). Using design motifs for inhibitory and excitatory regulations, three different oscillator designs were constructed and characterized: a two-switch negative-feedback oscillator, loosely analogous to the p53–Mdm2-feedback loop (Bar-Or et al, 2000); the same oscillator augmented with a positive-feedback loop, loosely analogous to a synthetic relaxation oscillator (Atkinson et al, 2003); and a three-switch ring oscillator analogous to the repressilator (Elowitz and Leibler, 2000).
DNA and RNA hybridization reactions (Figure 1B) can be assembled to create either an inhibitable switch (Figure 1A, right and bottom) with a threshold set by the total concentration of its DNA activator strand (Figure 1C, bottom), or an activatable switch (Figure 1A, left and top) with a threshold set by its DNA inhibitor strand concentration (Figure 1C, top). This threshold mechanism is analogous to biological threshold mechanisms such as ‘inhibitor ultrasensitivity' (Ferrell, 1996) and ‘molecular titration' (Buchler and Louis, 2008). Using these design motifs, we constructed a two-switch negative-feedback oscillator (Figure 1A, inset): RNA activator rA1 activates the production of RNA inhibitor rI2 by modulating switch Sw21, while RNA inhibitor rI2, in turn, inhibits the production of RNA activator rA1 by modulating switch Sw12. A total of seven DNA strands are used, in addition to the two enzymes, bacteriophage T7 RNA polymerase and Escherichia coli ribonuclease H. The fact that such a negative-feedback loop can lead to temporal oscillations can be seen from a mathematical model of transcriptional networks. Experimental results showed qualitative agreement with predicted oscillator behavior from simple model simulations.
The fully optimized system revealed five complete oscillation cycles with a nearly 50% amplitude swing (Figure 3A) until, after ∼20 h, the production rate could no longer be sustained in the batch reaction. Gel measurements verified oscillations in RNA concentrations and switch states (Figure 3B and C). However to our surprise, rather than oscillations with constant amplitude and constant mean, the RNA inhibitor concentration builds up after each cycle. An extended mathematical model that incorporated an interference reaction from ‘waste' product (Figure 3B and C) could qualitatively capture this behavior.
Using a new autoregulatory switch Sw11, we added a positive-feedback loop to the two-node oscillator to make an amplified negative feedback oscillator (Design II, Figure 1D). Further, we replaced the excitatory connection of Sw21 by a chain of two inhibitory connections, Sw23 and Sw31, to construct a three-switch ring oscillator (Design III, Figure 1D). All three oscillator designs could be tuned to reach the oscillatory regime in parameter space.
Reassuringly, our in vitro oscillators exhibit several design principles previously observed in vivo. (1) Introducing delay in a simple negative-feedback loop can help achieve stable oscillation (Novák and Tyson, 2008; Stricker et al, 2008). (2) The addition of a positive-feedback self-loop to a negative-feedback oscillator provides access to rich dynamics and improved tunability (Tsai et al, 2008). (3) Oscillations in biochemical ring oscillators (such as the repressilator) are sensitive to parameter asymmetry among individual components (Tuttle et al, 2005). (4) The saturation of degradation machinery and the management of waste products could play an important role.
However, several significant difficulties remain for predictive and robust oscillator performances: limited lifetime of closed batch reactions, interference from waste products, and asymmetry of switch components make quantitative modeling and predictio difficult. As a complementary approach to top-down view of systems biology, cell-free in vitro systems offer a valuable training ground to create and explore increasingly interesting and powerful information-based chemical systems (Simpson, 2006). In vitro oscillators could be used to orchestrate other chemical processes such as DNA nanomachines (Dittmer and Simmel, 2004) and to provide embedded controllers within prototype artificial cells (Noireaux and Libchaber, 2004; Griffiths and Tawfik, 2006).
The construction of synthetic biochemical circuits from simple components illuminates how complex behaviors can arise in chemistry and builds a foundation for future biological technologies. A simplified analog of genetic regulatory networks, in vitro transcriptional circuits, provides a modular platform for the systematic construction of arbitrary circuits and requires only two essential enzymes, bacteriophage T7 RNA polymerase and Escherichia coli ribonuclease H, to produce and degrade RNA signals. In this study, we design and experimentally demonstrate three transcriptional oscillators in vitro. First, a negative feedback oscillator comprising two switches, regulated by excitatory and inhibitory RNA signals, showed up to five complete cycles. To demonstrate modularity and to explore the design space further, a positive-feedback loop was added that modulates and extends the oscillatory regime. Finally, a three-switch ring oscillator was constructed and analyzed. Mathematical modeling guided the design process, identified experimental conditions likely to yield oscillations, and explained the system's robust response to interference by short degradation products. Synthetic transcriptional oscillators could prove valuable for systematic exploration of biochemical circuit design principles and for controlling nanoscale devices and orchestrating processes within artificial cells.
PMCID: PMC3063688  PMID: 21283141
cell free; in vitro; oscillation; synthetic biology; transcriptional circuits
18.  Dynamical Modeling of Collective Behavior from Pigeon Flight Data: Flock Cohesion and Dispersion 
PLoS Computational Biology  2012;8(3):e1002449.
Several models of flocking have been promoted based on simulations with qualitatively naturalistic behavior. In this paper we provide the first direct application of computational modeling methods to infer flocking behavior from experimental field data. We show that this approach is able to infer general rules for interaction, or lack of interaction, among members of a flock or, more generally, any community. Using experimental field measurements of homing pigeons in flight we demonstrate the existence of a basic distance dependent attraction/repulsion relationship and show that this rule is sufficient to explain collective behavior observed in nature. Positional data of individuals over time are used as input data to a computational algorithm capable of building complex nonlinear functions that can represent the system behavior. Topological nearest neighbor interactions are considered to characterize the components within this model. The efficacy of this method is demonstrated with simulated noisy data generated from the classical (two dimensional) Vicsek model. When applied to experimental data from homing pigeon flights we show that the more complex three dimensional models are capable of simulating trajectories, as well as exhibiting realistic collective dynamics. The simulations of the reconstructed models are used to extract properties of the collective behavior in pigeons, and how it is affected by changing the initial conditions of the system. Our results demonstrate that this approach may be applied to construct models capable of simulating trajectories and collective dynamics using experimental field measurements of herd movement. From these models, the behavior of the individual agents (animals) may be inferred.
Author Summary
The construction of mathematical models from experimental time-series data has been considered with some success in many areas of science and engineering, using the power of computer algorithms to build model structures and suitably tuning their parameters. When considering complex systems with nonlinear or collective behavior, computational models built from real data are the alternative to emulating the system as best as possible, since classic modeling approaches based on observation could prove difficult for complex dynamics. In this study, we provide a method to build models of collective dynamics from homing pigeon flight data. We show that our models follow the source dynamics well, and from them we are able to infer that significant collective behavior occurs in pigeon flights. Our results are consistent with the basic principles of previous hypotheses and models that have been proposed. Our approach serves as an initial outline towards the usage of experimental data to construct computational models to understand many complex phenomena with hypothesized collective behavior.
PMCID: PMC3315451  PMID: 22479176
19.  State reduction for network intervention in probabilistic Boolean networks 
Bioinformatics  2010;26(24):3098-3104.
Motivation: A key goal of studying biological systems is to design therapeutic intervention strategies. Probabilistic Boolean networks (PBNs) constitute a mathematical model which enables modeling, predicting and intervening in their long-run behavior using Markov chain theory. The long-run dynamics of a PBN, as represented by its steady-state distribution (SSD), can guide the design of effective intervention strategies for the modeled systems. A major obstacle for its application is the large state space of the underlying Markov chain, which poses a serious computational challenge. Hence, it is critical to reduce the model complexity of PBNs for practical applications.
Results: We propose a strategy to reduce the state space of the underlying Markov chain of a PBN based on a criterion that the reduction least distorts the proportional change of stationary masses for critical states, for instance, the network attractors. In comparison to previous reduction methods, we reduce the state space directly, without deleting genes. We then derive stationary control policies on the reduced network that can be naturally induced back to the original network. Computational experiments study the effects of the reduction on model complexity and the performance of designed control policies which is measured by the shift of stationary mass away from undesirable states, those associated with undesirable phenotypes. We consider randomly generated networks as well as a 17-gene gastrointestinal cancer network, which, if not reduced, has a 217 × 217 transition probability matrix. Such a dimension is too large for direct application of many previously proposed PBN intervention strategies.
Supplementary information: Supplementary information are available at Bioinformatics online.
PMCID: PMC3025721  PMID: 20956246
20.  Global Entrainment of Transcriptional Systems to Periodic Inputs 
PLoS Computational Biology  2010;6(4):e1000739.
This paper addresses the problem of providing mathematical conditions that allow one to ensure that biological networks, such as transcriptional systems, can be globally entrained to external periodic inputs. Despite appearing obvious at first, this is by no means a generic property of nonlinear dynamical systems. Through the use of contraction theory, a powerful tool from dynamical systems theory, it is shown that certain systems driven by external periodic signals have the property that all their solutions converge to a fixed limit cycle. General results are proved, and the properties are verified in the specific cases of models of transcriptional systems as well as constructs of interest in synthetic biology. A self-contained exposition of all needed results is given in the paper.
Author Summary
The activities of living organisms are governed by complex sets of biochemical reactions. Often, entrainment to certain external signals helps control the timing and sequencing of reactions. An important open problem is to understand the onset of entrainment and under what conditions it can be ensured in the presence of uncertainties, noise, and environmental variations. In this paper, we focus mainly on transcriptional systems, modeled by Ordinary Differential Equations. These are basic building blocks for more complex biochemical systems. However, the results that we obtain are of more generality. To illustrate this generality, and to emphasize the use of our techniques in synthetic biology, we discuss the entrainment of a Repressilator circuit and the synchronization of a network of Repressilators. We answer the following two questions: 1) What are the dynamical mechanisms that ensure the entrainment to periodic inputs in transcriptional modules? 2) Starting from natural systems, what properties can be used to design novel synthetic biological circuits that can be entrained? For some biological systems which are always “in contact” with a continuously changing environment, entrainment may be a “desired” property. Thus, answering the above two questions is of fundamental importance. While entrainment may appear obvious at first thought, it is not a generic property of nonlinear dynamical systems. The main result of our paper shows that, even if the transcriptional modules are modeled by nonlinear ODEs, they can be entrained by any (positive) periodic signal. Surprisingly, such a property is preserved if the system parameters are varied: entrainment is obtained independently of the particular biochemical conditions. We prove that combinations of the above transcriptional module also show the same property. Finally, we show how the developed tools can be applied to design synthetic biochemical systems guaranteed to exhibit entrainment.
PMCID: PMC2855316  PMID: 20418962
21.  Acoustic signatures of sound source-tract coupling 
Birdsong is a complex behavior, which results from the interaction between a nervous system and a biomechanical peripheral device. While much has been learned about how complex sounds are generated in the vocal organ, little has been learned about the signature on the vocalizations of the nonlinear effects introduced by the acoustic interactions between a sound source and the vocal tract. The variety of morphologies among bird species makes birdsong a most suitable model to study phenomena associated to the production of complex vocalizations. Inspired by the sound production mechanisms of songbirds, in this work we study a mathematical model of a vocal organ, in which a simple sound source interacts with a tract, leading to a delay differential equation. We explore the system numerically, and by taking it to the weakly nonlinear limit, we are able to examine its periodic solutions analytically. By these means we are able to explore the dynamics of oscillatory solutions of a sound source-tract coupled system, which are qualitatively different from those of a sound source-filter model of a vocal organ. Nonlinear features of the solutions are proposed as the underlying mechanisms of observed phenomena in birdsong, such as unilaterally produced “frequency jumps,” enhancement of resonances, and the shift of the fundamental frequency observed in heliox experiments.
PMCID: PMC3909991  PMID: 21599213
22.  Dynamics within the CD95 death-inducing signaling complex decide life and death of cells 
CD95-mediated apoptotic and NF-κB signaling were described by a simple kinetic model. We used a model reduction technique to reduce the number of reactions from 92 to 23 while maintaining a good model fit.p43-FLIP, which is generated at the CD95 DISC by procaspase-8 cleavage, was found to be the link between the CD95 DISC and the NF-κB pathway. P43-FLIP interacts with the IKK complex and leads to its activation.The CD95 DISC complex acts as a signal processor that diverges signals into the apoptotic and NF-κB pathways depending on the amounts of specific DISC proteins.Life/death decisions in CD95 signaling are determined by c-FLIPL and procaspase-8 in a non-linear way.
The CD95 protein (APO-1/Fas; Krammer et al, 2007) is a member of the death receptor family. Signal transduction of CD95 starts with the formation of the death-inducing signaling complex (DISC) detectable within seconds after receptor stimulation (Kischkel et al, 1995). The DISC consists of CD95, the adaptor molecule FADD, procaspase-8/10 and c-FLIPL/S/R (Muzio et al, 1996; Scaffidi et al, 1999; Sprick et al, 2002; Golks et al, 2005; Krammer et al, 2007). Procaspase-8 is converted at the DISC, in a series of autoproteolytic cleavage steps, to p43/p41 and p18, which leads to the activation of effector caspase-3 and demolition of the cell. Recently, experiments have demonstrated that CD95L also activates the induction of transcription factor NF-κB (Barnhart et al, 2004; Kreuz et al, 2004; Peter et al, 2007). It was shown that DED-containing proteins at the DISC, such as procaspase-8 and c-FLIP have a complex role in NF-κB activation (Chaudhary et al, 2000; Hu et al, 2000; Kreuz et al, 2004; Dohrman et al, 2005; Su et al, 2005). These findings motivated our systems biology approach and prompted us to determine whether CD95-mediated signaling should be considered a dynamic system, resulting in life/death decisions.
We observed simultaneous apoptosis and NF-κB induction on CD95 stimulation in HeLa cells stably overexpressing CD95–GFP (HeLa-CD95) using biochemical approaches and live-cell imaging. To understand the crosstalk between CD95-mediated apoptosis and NF-κB activation, we created a mathematical model of CD95 signaling. Our model assumes a trimerized ligand (L) that binds to a trimerized CD95 receptor (R) that can recruit three copies of FADD (F) leading to the DISC formation. Subsequently, DED-containing proteins, such as procaspase-8 (C8), c-FLIPL (FL) and c-FLIPS (FS) can bind to FADD. The order of protein binding gives rise to a combinatorial variety of intermediates, resulting either in the formation of the cleavage product of procaspase-8: p43/p41, or in the formation of the cleavage product of c-FLIPL: p43-FLIP. p43/p41 gives rise to signaling in the apoptotic branch of the model, whereas the cleavage product p43-FLIP triggers the activation of NF-κB. The model postulates that p43-FLIP interacts with the IKK complex leading to the phosphorylation of IκB (NF-κB·IκB·P), which entails its degradation and the translocation of p65 to the nucleus (NF-κB*). As a validation of the model topology, we confirmed experimentally that p43-FLIP interacts with the IKK complex and subsequently leads to its activation.
The complete model could be fitted well to a data set derived from quantitative western blots of a number of key proteins of the apoptotic and NF-κB pathways. However, we tested whether all the 92 reactions were required to reproduce the observed dynamics, as a small model would yield more reliable parameter estimates, which in turn would increase its usefulness as a predictive tool. To determine the most important interactions, we simplified the complete model in a step-wise manner obtaining a model of considerably lower complexity (Figure 5A, simplification steps are listed in Figure 5B). The final reduced model still approximated well the experimental data set (Figure 5D), whereas the number of reactions decreased from 92 to 23 (Figure 5C).
To better understand the interplay of DISC proteins in the determination of cell fate, we analyzed the activity of caspase-3 and NF-κB as a function of procaspase-8 and c-FLIPL levels (Figure 8A). We observed in our simulations that the decision over apoptosis and NF-κB is controlled by both proteins. Different scenarios occur that show combination or absence of either caspase-3 or NF-κB activity. The phase diagram shown in Figure 8A predicts that either increasing or decreasing the amount of c-FLIPL leads to a different signaling mode. We sought to validate this prediction by downregulating or overexpressing procaspase-8 and c-FLIPL, respectively, in HeLa-CD95 cells and measuring CD95-mediated signaling. In agreement with the phase diagram (Figure 8A), we observed that c-FLIPL overexpression resulted in a strong reduction of apoptosis (Figure 8D). Furthermore, we could further confirm by western blot analysis that the stable knockdown of c-FLIPL and procaspase-8 led to a reduction of the levels of p43-FLIP and phosphorylated IκBα after receptor stimulation (Figure 8C and D). In addition, to control the specificity of c-FLIP downregulation and further confirm the requirement of cleavage of c-FLIPL to p43-FLIP, we performed a reconstitution experiment in HeLa-CD95–c-FLIP-deficient cells (Figure 8E). Cells reconstituted with WT c-FLIPL were able to generate p43-FLIP and increased IκBα phosphorylation on CD95 stimulation. In contrast, cells reconstituted with the noncleavable mutant of c-FLIPL (D376E) did not show processing to p43-FLIP (Figure 8E; Supplementary Figure S9). Noticeably, as postulated by the model, this resulted in a strong reduction of the levels of IκBα phosphorylation on CD95 stimulation. Hence, by perturbing the ratio of procaspase-8 to c-FLIPL at the DISC, we directed the induction of apoptosis and NF-κB activation as predicted by our model. Taken together, we found that the DISC protein levels determine cell fate in a nonlinear manner, highlighting the role of signal processing within the DISC.
In this study, we propose, to the best of our knowledge, the first integrated kinetic model of CD95-mediated apoptosis and NF-κB signaling. This was achieved by integrating mechanistic knowledge of DISC assembly and caspase activation with a simple scheme of NF-κB activation. We observed that c-FLIPL levels crucially determine the balance between apoptotic and NF-κB signaling by shaping the dynamics of DISC assembly. Although this finding is based on experiments performed in cell lines, we expect that the nonlinear dynamics of DISC assembly is a generic systems property of life/death decision making in CD95 signaling pathways. This is especially important for understanding the regulation of cell death in physiologically relevant cells, such as cancer cells often showing resistance against death receptor-induced apoptosis.
This study explores the dilemma in cellular signaling that triggering of CD95 (Fas/APO-1) in some situations results in cell death and in others leads to the activation of NF-κB. We established an integrated kinetic mathematical model for CD95-mediated apoptotic and NF-κB signaling. Systematic model reduction resulted in a surprisingly simple model well approximating experimentally observed dynamics. The model postulates a new link between c-FLIPL cleavage in the death-inducing signaling complex (DISC) and the NF-κB pathway. We validated experimentally that CD95 stimulation resulted in an interaction of p43-FLIP with the IKK complex followed by its activation. Furthermore, we showed that the apoptotic and NF-κB pathways diverge already at the DISC. Model and experimental analysis of DISC formation showed that a subtle balance of c-FLIPL and procaspase-8 determines life/death decisions in a nonlinear manner. We present an integrated model describing the complex dynamics of CD95-mediated apoptosis and NF-κB signaling.
PMCID: PMC2858442  PMID: 20212524
apoptosis; CD95 signaling; DISC; model reduction; NF-κB
23.  Strategies for the Physiome Project 
Annals of biomedical engineering  2000;28(8):1043-1058.
The physiome is the quantitative description of the functioning organism in normal and pathophysiological states. The human physiome can be regarded as the virtual human. It is built upon the morphome, the quantitative description of anatomical structure, chemical and biochemical composition, and material properties of an intact organism, including its genome, proteome, cell, tissue, and organ structures up to those of the whole intact being. The Physiome Project is a multicentric integrated program to design, develop, implement, test and document, archive and disseminate quantitative information, and integrative models of the functional behavior of molecules, organelles, cells, tissues, organs, and intact organisms from bacteria to man. A fundamental and major feature of the project is the databasing of experimental observations for retrieval and evaluation. Technologies allowing many groups to work together are being rapidly developed. Internet II will facilitate this immensely. When problems are huge and complex, a particular working group can be expert in only a small part of the overall project. The strategies to be worked out must therefore include how to pull models composed of many submodules together even when the expertise in each is scattered amongst diverse institutions. The technologies of bioinformatics will contribute greatly to this effort. Developing and implementing code for large-scale systems has many problems. Most of the submodules are complex, requiring consideration of spatial and temporal events and processes. Submodules have to be linked to one another in a way that preserves mass balance and gives an accurate representation of variables in nonlinear complex biochemical networks with many signaling and controlling pathways. Microcompartmentalization vitiates the use of simplified model structures. The stiffness of the systems of equations is computationally costly. Faster computation is needed when using models as thinking tools and for iterative data analysis. Perhaps the most serious problem is the current lack of definitive information on kinetics and dynamics of systems, due in part to the almost total lack of databased observations, but also because, though we are nearly drowning in new information being published each day, either the information required for the modeling cannot be found or has never been obtained. “Simple” things like tissue composition, material properties, and mechanical behavior of cells and tissues are not generally available. The development of comprehensive models of biological systems is a key to pharmaceutics and drug design, for the models will become gradually better predictors of the results of interventions, both genomic and pharmaceutic. Good models will be useful in predicting the side effects and long term effects of drugs and toxins, and when the models are really good, to predict where genomic intervention will be effective and where the multiple redundancies in our biological systems will render a proposed intervention useless. The Physiome Project will provide the integrating scientific basis for the Genes to Health initiative, and make physiological genomics a reality applicable to whole organisms, from bacteria to man.
PMCID: PMC3425440  PMID: 11144666
Simulation analysis; Biological systems modeling; Complexity; Chaos; Dynamic Systems; Milieu interieure; Homeodynamics; Homeostasis; Metabolism; Control
24.  A Scalable Computational Framework for Establishing Long-Term Behavior of Stochastic Reaction Networks 
PLoS Computational Biology  2014;10(6):e1003669.
Reaction networks are systems in which the populations of a finite number of species evolve through predefined interactions. Such networks are found as modeling tools in many biological disciplines such as biochemistry, ecology, epidemiology, immunology, systems biology and synthetic biology. It is now well-established that, for small population sizes, stochastic models for biochemical reaction networks are necessary to capture randomness in the interactions. The tools for analyzing such models, however, still lag far behind their deterministic counterparts. In this paper, we bridge this gap by developing a constructive framework for examining the long-term behavior and stability properties of the reaction dynamics in a stochastic setting. In particular, we address the problems of determining ergodicity of the reaction dynamics, which is analogous to having a globally attracting fixed point for deterministic dynamics. We also examine when the statistical moments of the underlying process remain bounded with time and when they converge to their steady state values. The framework we develop relies on a blend of ideas from probability theory, linear algebra and optimization theory. We demonstrate that the stability properties of a wide class of biological networks can be assessed from our sufficient theoretical conditions that can be recast as efficient and scalable linear programs, well-known for their tractability. It is notably shown that the computational complexity is often linear in the number of species. We illustrate the validity, the efficiency and the wide applicability of our results on several reaction networks arising in biochemistry, systems biology, epidemiology and ecology. The biological implications of the results as well as an example of a non-ergodic biological network are also discussed.
Author Summary
In many biological disciplines, computational modeling of interaction networks is the key for understanding biological phenomena. Such networks are traditionally studied using deterministic models. However, it has been recently recognized that when the populations are small in size, the inherent random effects become significant and to incorporate them, a stochastic modeling paradigm is necessary. Hence, stochastic models of reaction networks have been broadly adopted and extensively used. Such models, for instance, form a cornerstone for studying heterogeneity in clonal cell populations. In biological applications, one is often interested in knowing the long-term behavior and stability properties of reaction networks even with incomplete knowledge of the model parameters. However for stochastic models, no analytical tools are known for this purpose, forcing many researchers to use a simulation-based approach, which is highly unsatisfactory. To address this issue, we develop a theoretical and computational framework for determining the long-term behavior and stability properties for stochastic reaction networks. Our approach is based on a mixture of ideas from probability theory, linear algebra and optimization theory. We illustrate the broad applicability of our results by considering examples from various biological areas. The biological implications of our results are discussed as well.
PMCID: PMC4072526  PMID: 24968191
25.  A new flexible plug and play scheme for modeling, simulating, and predicting gastric emptying 
In-silico models that attempt to capture and describe the physiological behavior of biological organisms, including humans, are intrinsically complex and time consuming to build and simulate in a computing environment. The level of detail of description incorporated in the model depends on the knowledge of the system’s behavior at that level. This knowledge is gathered from the literature and/or improved by knowledge obtained from new experiments. Thus model development is an iterative developmental procedure. The objective of this paper is to describe a new plug and play scheme that offers increased flexibility and ease-of-use for modeling and simulating physiological behavior of biological organisms.
This scheme requires the modeler (user) first to supply the structure of the interacting components and experimental data in a tabular format. The behavior of the components described in a mathematical form, also provided by the modeler, is externally linked during simulation. The advantage of the plug and play scheme for modeling is that it requires less programming effort and can be quickly adapted to newer modeling requirements while also paving the way for dynamic model building.
As an illustration, the paper models the dynamics of gastric emptying behavior experienced by humans. The flexibility to adapt the model to predict the gastric emptying behavior under varying types of nutrient infusion in the intestine (ileum) is demonstrated. The predictions were verified with a human intervention study. The error in predicting the half emptying time was found to be less than 6%.
A new plug-and-play scheme for biological systems modeling was developed that allows changes to the modeled structure and behavior with reduced programming effort, by abstracting the biological system into a network of smaller sub-systems with independent behavior. In the new scheme, the modeling and simulation becomes an automatic machine readable and executable task.
PMCID: PMC4080776  PMID: 24917054
Modeling; Gastric emptying; Functional modules; Feedback loop

Results 1-25 (1016549)