Search

We demonstrate a modeling and computational framework that allows for rapid screening of thousands of potential network designs for particular dynamic behavior. To illustrate this capability we consider the problem of hysteresis, a prerequisite for construction of robust bistable switches and hence a cornerstone for construction of more complex synthetic circuits. We evaluate and rank most three node networks according to their ability to robustly exhibit hysteresis where robustness is measured with respect to parameters over multiple dynamic phenotypes. Focusing on the highest ranked networks, we demonstrate how additional robustness and design constraints can be applied. We compare our results to more traditional methods based on specific parameterization of ordinary differential equation models and demonstrate a strong qualitative match at a small fraction of the computational cost., Author summary A major challenge in the domains of systems and synthetic biology is an inability to efficiently predict function(s) of complex networks. This work demonstrates a modeling and computational framework that allows for a mathematically justifiable rigorous screening of thousands of potential network designs for a wide variety of dynamical behavior. We screen all 3-node genetic networks and rank them based on their ability to act as an inducible bistable switch. Our results are summarized in a searchable database that can be used to construct robust switches. The ability to quickly screen thousands of designs significantly reduces the set of viable designs and allows synthetic biologists to focus their experimental and more traditional modeling tools to this much smaller set.

The Hippo-YAP/TAZ signaling pathway plays a critical role in tissue homeostasis, tumorigenesis, and degeneration disorders. The regulation of YAP/TAZ levels is controlled by a complex regulatory network, where several feedback loops have been identified. However, it remains elusive how these feedback loops contain the YAP/TAZ levels and maintain the system in a healthy physiological state or trap the system in pathological conditions. Here, a mathematical model was developed to represent the YAP/TAZ regulatory network. Through theoretical analyses, three distinct states that designate the one physiological and two pathological outcomes were found. The transition from the physiological state to the two pathological states is mechanistically controlled by coupled bidirectional bistable switches, which are robust to parametric variation and stochastic fluctuations at the molecular level. This work provides a mechanistic understanding of the regulation and dysregulation of YAP/TAZ levels in tissue state transitions., Author summary Tissue development and homeostasis require well-controlled cell proliferation. Lack of this control could lead to degenerative or tumorigenic diseases. Signaling pathways have been explored in promoting or inhibiting these diseases. The Hippo signaling pathway is one of these, which has been found to control tissue homeostasis and organ size through cell proliferation and apoptosis, as evidenced by extensive experimental data. However, the question remains of how tissue can transition from a homeostatic state to either a degenerative or tumorigenic state. By theoretically analyzing a mathematical model of its regulatory network, we present a mechanism that underlies Hippo signaling to control tissue transition from a homeostatic state to a disease state. This provides us with a mechanistic understanding of how the parts of the regulatory network are coordinated for the transitions between the homeostasis state and the disease states. In addition, we looked at the role of system noise and found that it could promote the transition to one of the disease states. Our model allows for experimental hypotheses to be generated and could lead to the development of therapeutic strategies by targeting the Hippo signaling pathway.

Once biological systems are modeled by regulatory networks, the next step is to include external stimuli, which model the experimental possibilities to affect the activity level of certain network’s nodes, in a mathematical framework. Then, this framework can be interpreted as a mathematical optimal control framework such that optimization algorithms can be used to determine external stimuli which cause a desired switch from an initial state of the network to another final state. These external stimuli are the intervention points for the corresponding biological experiment to obtain the desired outcome of the considered experiment. In this work, the model of regulatory networks is extended to controlled regulatory networks. For this purpose, external stimuli are considered which can affect the activity of the network’s nodes by activation or inhibition. A method is presented how to calculate a selection of external stimuli which causes a switch between two different steady states of a regulatory network. A software solution based on Jimena and Mathworks Matlab is provided. Furthermore, numerical examples are presented to demonstrate application and scope of the software on networks of 4 nodes, 11 nodes and 36 nodes. Moreover, we analyze the aggregation of platelets and the behavior of a basic T-helper cell protein-protein interaction network and its maturation towards Th0, Th1, Th2, Th17 and Treg cells in accordance with experimental data., Author summary Organisms can be seen as molecular networks being able to react on external stimuli. Experiments are performed to understand the underlying regulating mechanisms within the molecular network. A common purpose for these efforts is to reveal mechanisms with which the molecular networks can be affected to achieve a desired behavior. To cover the complexity of life these models of molecular networks often need to be quite huge and need to have many cross connections between the different agents of the network that regulate the behavior of the network. A useful tool to structure this complexity are mathematical methods. Once the model based on experiments is set up the experimental data can be further processed by mathematical methods. As experiments are cumbersome, the present work provides a framework that can be used to systematically figure out intervention points in molecular networks to cause a desired effect. In this way promising intervention strategies can be obtained. For instance the process of obtaining new drugs for pharmacological modulation can be shortened as in the best case the calculated intervention strategy just has to be validated with one experiment and the time consuming procedure of searching an intervention strategy with several experiments can be saved.

Control of axial polarity during regeneration is a crucial open question. We developed a quantitative model of regenerating planaria, which elucidates self-assembly mechanisms of morphogen gradients required for robust body-plan control. The computational model has been developed to predict the fraction of heteromorphoses expected in a population of regenerating planaria fragments subjected to different treatments, and for fragments originating from different regions along the anterior-posterior and medio-lateral axis. This allows for a direct comparison between computational and experimental regeneration outcomes. Vector transport of morphogens was identified as a fundamental requirement to account for virtually scale-free self-assembly of the morphogen gradients observed in planarian homeostasis and regeneration. The model correctly describes altered body-plans following many known experimental manipulations, and accurately predicts outcomes of novel cutting scenarios, which we tested. We show that the vector transport field coincides with the alignment of nerve axons distributed throughout the planarian tissue, and demonstrate that the head-tail axis is controlled by the net polarity of neurons in a regenerating fragment. This model provides a comprehensive framework for mechanistically understanding fundamental aspects of body-plan regulation, and sheds new light on the role of the nervous system in directing growth and form., Author summary Understanding how large-scale anatomy emerges from the activity of cellular pathways is a key goal of evolutionary developmental biology. Elucidating the rules of body-wide morphogenesis is especially essential for transitioning molecular signaling data at the cellular level into advances in regenerative biomedicine. We constructed and analyzed a comprehensive, multiscale computational model to explain the determination of axial polarity during planarian regeneration. Uniquely, our model explains the various head-tail patterning outcomes of a wide range of molecular and physiological manipulations. Testing the novel predictions of this model revealed the nervous system as an instructive regulator of axial patterning.

We present a new modeling and computational tool that computes rigorous summaries of network dynamics over large sets of parameter values. These summaries, organized in a database, can be searched for observed dynamics, e.g., bistability and hysteresis, to discover parameter regimes over which they are supported. We illustrate our approach on several networks underlying the restriction point of the cell cycle in humans and yeast. We rank networks by how robustly they support hysteresis, which is the observed phenotype. We find that the best 6-node human network and the yeast network share similar topology and robustness of hysteresis, in spite of having no homology between the corresponding nodes of the network. Our approach provides a new tool linking network structure and dynamics., Author summary To summarize our understanding of how genes, their products and other cellular actors interact with each other, we often employ networks to describe their interactions. However, networks do not fully specify how the underlying biological system behaves in different conditions, nor how such response evolves in time. We present a new modeling and computational approach that allows us to compute and collect summaries of network dynamics for large sets of parameter values. We can then search these summaries for all observed behavior. We illustrate our approach on networks that govern entry to the cell cycle in humans and yeast. We rank networks based on how robustly they exhibit the experimentally observed behavior of hysteresis. We find similarities in network structure of the best ranked networks in yeast and humans, which are not explained by a common ancestry. Our approach provides a tool linking network structure and the behavior of the underlying system.

Current state-of-the-art mathematical models to investigate complex biological processes, in particular liver-associated pathologies, have limited expansiveness, flexibility, representation of integrated regulation and rely on the availability of detailed kinetic data. We generated the SteatoNet, a multi-pathway, multi-tissue model and in silico platform to investigate hepatic metabolism and its associated deregulations. SteatoNet is based on object-oriented modelling, an approach most commonly applied in automotive and process industries, whereby individual objects correspond to functional entities. Objects were compiled to feature two novel hepatic modelling aspects: the interaction of hepatic metabolic pathways with extra-hepatic tissues and the inclusion of transcriptional and post-transcriptional regulation. SteatoNet identification at normalised steady state circumvents the need for constraining kinetic parameters. Validation and identification of flux disturbances that have been proven experimentally in liver patients and animal models highlights the ability of SteatoNet to effectively describe biological behaviour. SteatoNet identifies crucial pathway branches (transport of glucose, lipids and ketone bodies) where changes in flux distribution drive the healthy liver towards hepatic steatosis, the primary stage of non-alcoholic fatty liver disease. Cholesterol metabolism and its transcription regulators are highlighted as novel steatosis factors. SteatoNet thus serves as an intuitive in silico platform to identify systemic changes associated with complex hepatic metabolic disorders., Author Summary In this article we present SteatoNet, the most comprehensive computational network of hepatic metabolism and the interaction of the liver with extra-hepatic tissues. Generation of the SteatoNet involved the application of engineering strategies to resolve prevalent drawbacks in biological modelling and thus effectively understand basic biology. SteatoNet does not require detailed kinetic parameters, behaves representatively of biological observations and portrays systemic interactions. SteatoNet is simple and flexible, which is an important advantage over current computational models of complex metabolic disorders associated with the liver. To demonstrate the utility of SteatoNet as a hypotheses-generation tool, we studied candidate mechanisms of hepatic fat accumulation, the key characteristic of non-alcoholic fatty liver disease (NAFLD), a common but poorly understood metabolic disorder influenced by genetic and lifestyle factors. Our data describe NAFLD as a network disease, with deregulated pathways such as inter-tissue transport of glucose, lipid and ketone bodies, cholesterol metabolism and “stiff” focal points in the fat tissue. These results indicate the involvement of systemic metabolic deregulations in the transformation of healthy to fatty liver in NAFLD. The SteatoNet highlights the utility of engineering approaches in systems biology to aid in solving complex biological questions.

Many developmental, physiological, and behavioral processes depend on the precise expression of genes in space and time. Such spatiotemporal gene expression phenotypes arise from the binding of sequence-specific transcription factors (TFs) to DNA, and from the regulation of nearby genes that such binding causes. These nearby genes may themselves encode TFs, giving rise to a transcription factor network (TFN), wherein nodes represent TFs and directed edges denote regulatory interactions between TFs. Computational studies have linked several topological properties of TFNs — such as their degree distribution — with the robustness of a TFN's gene expression phenotype to genetic and environmental perturbation. Another important topological property is assortativity, which measures the tendency of nodes with similar numbers of edges to connect. In directed networks, assortativity comprises four distinct components that collectively form an assortativity signature. We know very little about how a TFN's assortativity signature affects the robustness of its gene expression phenotype to perturbation. While recent theoretical results suggest that increasing one specific component of a TFN's assortativity signature leads to increased phenotypic robustness, the biological context of this finding is currently limited because the assortativity signatures of real-world TFNs have not been characterized. It is therefore unclear whether these earlier theoretical findings are biologically relevant. Moreover, it is not known how the other three components of the assortativity signature contribute to the phenotypic robustness of TFNs. Here, we use publicly available DNaseI-seq data to measure the assortativity signatures of genome-wide TFNs in 41 distinct human cell and tissue types. We find that all TFNs share a common assortativity signature and that this signature confers phenotypic robustness to model TFNs. Lastly, we determine the extent to which each of the four components of the assortativity signature contributes to this robustness., Author Summary The cells of living organisms do not concurrently express their entire complement of genes. Instead, they regulate their gene expression, and one consequence of this is the potential for different cells to adopt different stable gene expression patterns. For example, the development of an embryo necessitates that cells alter their gene expression patterns in order to differentiate. These gene expression phenotypes are largely robust to genetic mutation, and one source of this robustness may reside in the network structure of interacting molecules that underlie genetic regulation. Theoretical studies of regulatory networks have linked network structure to robustness; however, it is also necessary to more extensively characterize real-world regulatory networks in order to understand which structural properties may be biologically meaningful. We recently used theoretical models to show that a particular structural property, degree assortativity, is linked to robustness. Here, we measure the assortativity of human regulatory networks in 41 distinct cell and tissue types. We then develop a theoretical framework to explore how this structural property affects robustness, and we find that the gene expression phenotypes of human regulatory networks are more robust than expected by chance alone.

Oligomers of length k, or k-mers, are convenient and widely used features for modeling the properties and functions of DNA and protein sequences. However, k-mers suffer from the inherent limitation that if the parameter k is increased to resolve longer features, the probability of observing any specific k-mer becomes very small, and k-mer counts approach a binary variable, with most k-mers absent and a few present once. Thus, any statistical learning approach using k-mers as features becomes susceptible to noisy training set k-mer frequencies once k becomes large. To address this problem, we introduce alternative feature sets using gapped k-mers, a new classifier, gkm-SVM, and a general method for robust estimation of k-mer frequencies. To make the method applicable to large-scale genome wide applications, we develop an efficient tree data structure for computing the kernel matrix. We show that compared to our original kmer-SVM and alternative approaches, our gkm-SVM predicts functional genomic regulatory elements and tissue specific enhancers with significantly improved accuracy, increasing the precision by up to a factor of two. We then show that gkm-SVM consistently outperforms kmer-SVM on human ENCODE ChIP-seq datasets, and further demonstrate the general utility of our method using a Naïve-Bayes classifier. Although developed for regulatory sequence analysis, these methods can be applied to any sequence classification problem., Author Summary Genomic regulatory elements (enhancers, promoters, and insulators) control the expression of their target genes and are widely believed to play a key role in human development and disease by altering protein concentrations. A fundamental step in understanding enhancers is the development of DNA sequence-based models to predict the tissue specific activity of regulatory elements. Such models facilitate both the identification of the molecular pathways which impinge on enhancer activity through direct transcription factor binding, and the direct evaluation of the impact of specific common or rare genetic variants on enhancer function. We have previously developed a successful sequence-based model for enhancer prediction using a k-mer support vector machine (kmer-SVM). Here, we address a significant limitation of the kmer-SVM approach and present an alternative method using gapped k-mers (gkm-SVM) which exhibits dramatically improved accuracy in all test cases. While we focus on enhancers and transcription factor binding, our method can be applied to improve a much broader class of sequence analysis problems, including proteins and RNA. In addition, we expect that most k-mer based methods can be significantly improved by simply using the generalized k-mer count method that we present in this paper. We believe this improved model will enable significant contributions to our understanding of the human regulatory system.

Identifying master regulators of biological processes and mapping their downstream gene networks are key challenges in systems biology. We developed a computational method, called iRegulon, to reverse-engineer the transcriptional regulatory network underlying a co-expressed gene set using cis-regulatory sequence analysis. iRegulon implements a genome-wide ranking-and-recovery approach to detect enriched transcription factor motifs and their optimal sets of direct targets. We increase the accuracy of network inference by using very large motif collections of up to ten thousand position weight matrices collected from various species, and linking these to candidate human TFs via a motif2TF procedure. We validate iRegulon on gene sets derived from ENCODE ChIP-seq data with increasing levels of noise, and we compare iRegulon with existing motif discovery methods. Next, we use iRegulon on more challenging types of gene lists, including microRNA target sets, protein-protein interaction networks, and genetic perturbation data. In particular, we over-activate p53 in breast cancer cells, followed by RNA-seq and ChIP-seq, and could identify an extensive up-regulated network controlled directly by p53. Similarly we map a repressive network with no indication of direct p53 regulation but rather an indirect effect via E2F and NFY. Finally, we generalize our computational framework to include regulatory tracks such as ChIP-seq data and show how motif and track discovery can be combined to map functional regulatory interactions among co-expressed genes. iRegulon is available as a Cytoscape plugin from http://iregulon.aertslab.org., Author Summary Gene regulatory networks control developmental, homeostatic, and disease processes by governing precise levels and spatio-temporal patterns of gene expression. Determining their topology can provide mechanistic insight into these processes. Gene regulatory networks consist of interactions between transcription factors and their direct target genes. Each regulatory interaction represents the binding of the transcription factor to a specific DNA binding site near its target gene. Here we present a computational method, called iRegulon, to identify master regulators and direct target genes in a human gene signature, i.e. a set of co-expressed genes. iRegulon relies on the analysis of the regulatory sequences around each gene in the gene set to detect enriched TF motifs or ChIP-seq peaks, using databases of nearly 10.000 TF motifs and 1000 ChIP-seq data sets or “tracks”. Next, it associates enriched motifs and tracks with candidate transcription factors and determines the optimal subset of direct target genes. We validate iRegulon on ENCODE data, and use it in combination with RNA-seq and ChIP-seq data to map a p53 downstream network with new predicted co-factors and targets. iRegulon is available as a Cytoscape plugin, supporting human, mouse, and Drosophila genes, and provides access to hundreds of cancer-related TF-target subnetworks or “regulons”.

Gene regulatory networks are a crucial aspect of systems biology in describing molecular mechanisms of the cell. Various computational models rely on random gene selection to infer such networks from microarray data. While incorporation of prior knowledge into data analysis has been deemed important, in practice, it has generally been limited to referencing genes in probe sets and using curated knowledge bases. We investigate the impact of augmenting microarray data with semantic relations automatically extracted from the literature, with the view that relations encoding gene/protein interactions eliminate the need for random selection of components in non-exhaustive approaches, producing a more accurate model of cellular behavior. A genetic algorithm is then used to optimize the strength of interactions using microarray data and an artificial neural network fitness function. The result is a directed and weighted network providing the individual contribution of each gene to its target. For testing, we used invasive ductile carcinoma of the breast to query the literature and a microarray set containing gene expression changes in these cells over several time points. Our model demonstrates significantly better fitness than the state-of-the-art model, which relies on an initial random selection of genes. Comparison to the component pathways of the KEGG Pathways in Cancer map reveals that the resulting networks contain both known and novel relationships. The p53 pathway results were manually validated in the literature. 60% of non-KEGG relationships were supported (74% for highly weighted interactions). The method was then applied to yeast data and our model again outperformed the comparison model. Our results demonstrate the advantage of combining gene interactions extracted from the literature in the form of semantic relations with microarray analysis in generating contribution-weighted gene regulatory networks. This methodology can make a significant contribution to understanding the complex interactions involved in cellular behavior and molecular physiology., Author Summary We have developed a methodology that combines standard computational analysis of gene expression data with knowledge in the literature to identify pathways of gene and protein interactions. We extract the knowledge from PubMed citations using a tool (SemRep) that identifies specific relationships between genes or proteins. We string together networks of individual interactions that are found within citations that refer to the target pathways. Upon this skeleton of interactions, we calculate the weight of the interaction with the gene expression data captured over multiple time points using state-of-the-art analysis algorithms. Not surprisingly, this approach of combining prior knowledge into the analysis process significantly improves the performance of the analysis. This work is most significant as an example of how the wealth of textual data related to gene interactions can be incorporated into computational analysis, not solely to identify this type of pathway (a gene regulatory network) but for any type of similar biological problem.

The TGF-β/Smad signaling system decreases its activity through strong negative regulation. Several molecular mechanisms of negative regulation have been published, but the relative impact of each mechanism on the overall system is unknown. In this work, we used computational and experimental methods to assess multiple negative regulatory effects on Smad signaling in HaCaT cells. Previously reported negative regulatory effects were classified by time-scale: degradation of phosphorylated R-Smad and I-Smad-induced receptor degradation were slow-mode effects, and dephosphorylation of R-Smad was a fast-mode effect. We modeled combinations of these effects, but found no combination capable of explaining the observed dynamics of TGF-β/Smad signaling. We then proposed a negative feedback loop with upregulation of the phosphatase PPM1A. The resulting model was able to explain the dynamics of Smad signaling, under both short and long exposures to TGF-β. Consistent with this model, immuno-blots showed PPM1A levels to be significantly increased within 30 min after TGF-β stimulation. Lastly, our model was able to resolve an apparent contradiction in the published literature, concerning the dynamics of phosphorylated R-Smad degradation. We conclude that the dynamics of Smad negative regulation cannot be explained by the negative regulatory effects that had previously been modeled, and we provide evidence for a new negative feedback loop through PPM1A upregulation. This work shows that tight coupling of computational and experiments approaches can yield improved understanding of complex pathways., Singapore-MIT Alliance, Mechanobiology Institute, Singapore, Institute of Bioengineering and Nanotechnology (Singapore), National University of Singapore (NUS Graduate School for Integrative Sciences and Engineering scholar), Singapore-MIT Alliance for Research and Technology

The Drosophila eggshell constitutes a remarkable system for the study of epithelial patterning, both experimentally and through computational modeling. Dorsal eggshell appendages arise from specific regions in the anterior follicular epithelium that covers the oocyte: two groups of cells expressing broad (roof cells) bordered by rhomboid expressing cells (floor cells). Despite the large number of genes known to participate in defining these domains and the important modeling efforts put into this developmental system, key patterning events still lack a proper mechanistic understanding and/or genetic basis, and the literature appears to conflict on some crucial points. We tackle these issues with an original, discrete framework that considers single-cell models that are integrated to construct epithelial models. We first build a phenomenological model that reproduces wild type follicular epithelial patterns, confirming EGF and BMP signaling input as sufficient to establish the major features of this patterning system within the anterior domain. Importantly, this simple model predicts an instructive juxtacrine signal linking the roof and floor domains. To explore this prediction, we define a mechanistic model that integrates the combined effects of cellular genetic networks, cell communication and network adjustment through developmental events. Moreover, we focus on the anterior competence region, and postulate that early BMP signaling participates with early EGF signaling in its specification. This model accurately simulates wild type pattern formation and is able to reproduce, with unprecedented level of precision and completeness, various published gain-of-function and loss-of-function experiments, including perturbations of the BMP pathway previously seen as conflicting results. The result is a coherent model built upon rules that may be generalized to other epithelia and developmental systems., Author Summary We propose a modeling approach to Drosophila egg appendage development. Relying on a thorough compilation of published data, our model comprises cellular genetic networks and cell-cell communication. It proves to be extremely robust by faithfully replicating the wild type pattern as well as patterns arising from mutational perturbations. Our model proposes a new hypothesis for the definition of the anterior competence region, through early posterior EGF activity in conjunction with an early BMP signaling event, which may reconcile apparently conflicting experimental results. Furthermore, this is the first model that explicitly includes the removal of Grk, possibly through vitelline membrane formation, as a central event in patterning the follicular epithelium. Importantly, the model predicts the requirement of an instructive juxtacrine signal to specify the different regions of the developing appendages. Finally, our modeling framework may be applied to other developing epithelia and biological systems.

Biological protein interactions networks such as signal transduction or gene transcription networks are often treated as modular, allowing motifs to be analyzed in isolation from the rest of the network. Modularity is also a key assumption in synthetic biology, where it is similarly expected that when network motifs are combined together, they do not lose their essential characteristics. However, the interactions that a network module has with downstream elements change the dynamical equations describing the upstream module and thus may change the dynamic and static properties of the upstream circuit even without explicit feedback. In this work we analyze the behavior of a ubiquitous motif in gene transcription and signal transduction circuits: the switch. We show that adding an additional downstream component to the simple genetic toggle switch changes its dynamical properties by changing the underlying potential energy landscape, and skewing it in favor of the unloaded side, and in some situations adding loads to the genetic switch can also abrogate bistable behavior. We find that an additional positive feedback motif found in naturally occurring toggle switches could tune the potential energy landscape in a desirable manner. We also analyze autocatalytic signal transduction switches and show that a ubiquitous positive feedback switch can lose its switch-like properties when connected to a downstream load. Our analysis underscores the necessity of incorporating the effects of downstream components when understanding the physics of biochemical network motifs, and raises the question as to how these effects are managed in real biological systems. This analysis is particularly important when scaling synthetic networks to more complex organisms., Author Summary Cells rely on complex networks of protein-protein interactions in order to carry out life functions. Scientists believe that these networks are organized in a modular fashion; that is they are made up of functionally distinct parts like an electronic circuit. Modularity implies that just as we put together electronic parts to make an amplifier that we can use in many different circuits, we can put together biochemical reactions to make an amplifier, or a switch or an oscillator, which perform the same function in different organisms. This assumption is important in synthetic biology, where we engineer and assemble synthetic genetic circuits in living organisms in a modular fashion. We show that for important modules like genetic and signaling switches, the assumption of modularity has a crucial limitation. We show that if one simply connects a biological switch to another downstream circuit, the presence of the connection changes the operation of the switch, which in some cases may stop behaving like a switch. Our work underscores the importance of taking into account these downstream connections and suggests that natural systems may be balancing some of these components in order to ensure that despite diversity, modules continue to work in different systems with fidelity.

Circadian clocks are biological timekeepers that allow living cells to time their activity in anticipation of predictable daily changes in light and other environmental factors. The complexity of the circadian clock in higher plants makes it difficult to understand the role of individual genes or molecular interactions, and mathematical modelling has been useful in guiding clock research in model organisms such as Arabidopsis thaliana. We present a model of the circadian clock in Arabidopsis, based on a large corpus of published time course data. It appears from experimental evidence in the literature that most interactions in the clock are repressive. Hence, we remove all transcriptional activation found in previous models of this system, and instead extend the system by including two new components, the morning-expressed activator RVE8 and the nightly repressor/activator NOX. Our modelling results demonstrate that the clock does not need a large number of activators in order to reproduce the observed gene expression patterns. For example, the sequential expression of the PRR genes does not require the genes to be connected as a series of activators. In the presented model, transcriptional activation is exclusively the task of RVE8. Predictions of how strongly RVE8 affects its targets are found to agree with earlier interpretations of the experimental data, but generally we find that the many negative feedbacks in the system should discourage intuitive interpretations of mutant phenotypes. The dynamics of the clock are difficult to predict without mathematical modelling, and the clock is better viewed as a tangled web than as a series of loops., Author Summary Like most living organisms, plants are dependent on sunlight, and evolution has endowed them with an internal clock by which they can predict sunrise and sunset. The clock consists of many genes that control each other in a complex network, leading to daily oscillations in protein levels. The interactions between genes can be positive or negative, causing target genes to be turned on or off. By constructing mathematical models that incorporate our knowledge of this network, we can interpret experimental data by comparing with results from the models. Any discrepancy between experimental data and model predictions will highlight where we are lacking in understanding. We compiled more than 800 sets of measured data from published articles about the clock in the model organism thale cress (Arabidopsis thaliana). Using these data, we constructed a mathematical model which compares favourably with previous models for simulating the clock. We used our model to investigate the role of positive interactions between genes, whether they are necessary for the function of the clock and if they can be identified in the model.

We explore a framework to model the dose response of allosteric multisite phosphorylation proteins using a single auxiliary variable. This reduction can closely replicate the steady state behavior of detailed multisite systems such as the Monod-Wyman-Changeux allosteric model or rule-based models. Optimal ultrasensitivity is obtained when the activation of an allosteric protein by its individual sites is concerted and redundant. The reduction makes this framework useful for modeling and analyzing biochemical systems in practical applications, where several multisite proteins may interact simultaneously. As an application we analyze a newly discovered checkpoint signaling pathway in budding yeast, which has been proposed to measure cell growth by monitoring signals generated at sites of plasma membrane growth. We show that the known components of this pathway can form a robust hysteretic switch. In particular, this system incorporates a signal proportional to bud growth or size, a mechanism to read the signal, and an all-or-none response triggered only when the signal reaches a threshold indicating that sufficient growth has occurred., Author Summary A large number of proteins in the cell are modified post-translationally by phosphorylation at multiple specific locations. This can help bring about interesting dynamical behaviors such as bistability or all-or-none responses to stimuli. Such behaviors are in turn important for cellular decision-making, differentiation, or the regulation of cellular processes. In this paper we propose the use of a specific technique for modeling allosteric multisite proteins, which can be thought of as the reduced version of a more detailed set of reactions. After validating this technique computationally by comparison to more detailed systems, we apply it to a new model of a signal transduction pathway with several multisite proteins. The model is concerned with a mechanism in budding yeast that is thought to measure the extent of daughter bud growth and send a signal to initiate mitosis when sufficient growth has occurred. Using the given framework we derive an analytically tractable model that creates the desired all-or-none signal. Overall, we give quantitative support to a newly discovered biochemical pathway, and we show the utility of the new modeling framework in the context of a realistic biological problem.

The mucosa of the intestinal tract represents a finely tuned system where tissue structure strongly influences, and is turn influenced by, its function as both an absorptive surface and a defensive barrier. Mucosal architecture and histology plays a key role in the diagnosis, characterization and pathophysiology of a host of gastrointestinal diseases. Inflammation is a significant factor in the pathogenesis in many gastrointestinal diseases, and is perhaps the most clinically significant control factor governing the maintenance of the mucosal architecture by morphogenic pathways. We propose that appropriate characterization of the role of inflammation as a controller of enteric mucosal tissue patterning requires understanding the underlying cellular and molecular dynamics that determine the epithelial crypt-villus architecture across a range of conditions from health to disease. Towards this end we have developed the Spatially Explicit General-purpose Model of Enteric Tissue (SEGMEnT) to dynamically represent existing knowledge of the behavior of enteric epithelial tissue as influenced by inflammation with the ability to generate a variety of pathophysiological processes within a common platform and from a common knowledge base. In addition to reproducing healthy ileal mucosal dynamics as well as a series of morphogen knock-out/inhibition experiments, SEGMEnT provides insight into a range of clinically relevant cellular-molecular mechanisms, such as a putative role for Phosphotase and tensin homolog/phosphoinositide 3-kinase (PTEN/PI3K) as a key point of crosstalk between inflammation and morphogenesis, the protective role of enterocyte sloughing in enteric ischemia-reperfusion and chronic low level inflammation as a driver for colonic metaplasia. These results suggest that SEGMEnT can serve as an integrating platform for the study of inflammation in gastrointestinal disease., Author Summary Mucosal histology plays a key role in the diagnosis, characterization and propagation of a host of gastrointestinal diseases, and the development of computational models capable of producing spatial architecture comparable to histology will enhance the evaluation of hypotheses for those diseases. Inflammation is a significant factor in the pathogenesis of a series of gastrointestinal diseases, and affects the maintenance of the mucosal architecture by morphogenic pathways. We have developed the Spatially Explicit General-purpose Model of Enteric Tissue (SEGMEnT) to investigate the behavior of enteric epithelial tissue as influenced by inflammation. SEGMEnT integrates cellular and molecular pathways governing morphogenesis and inflammation to generate a variety of pathophysiological processes from a common platform and knowledge base. Beyond reproducing healthy and disease ileal mucosal dynamics, SEGMEnT provides insight into a range of clinically relevant cellular-molecular mechanisms, including a novel putative role for Phosphotase and tensin homolog/phosphoinositide 3-kinase (PTEN/PI3K) as a key point of crosstalk between enteric inflammation and morphogenesis, the protective role of enterocyte sloughing in enteric ischemia-reperfusion and the mechanism of chronic low level inflammation as a driver for colonic metaplasia. These results suggest that SEGMEnT can serve as an integrating platform for the study of inflammation in gastrointestinal disease.

Cells generally adapt to environmental changes by first exhibiting an immediate response and then gradually returning to their original state to achieve homeostasis. Although simple network motifs consisting of a few genes have been shown to exhibit such adaptive dynamics, they do not reflect the complexity of real cells, where the expression of a large number of genes activates or represses other genes, permitting adaptive behaviors. Here, we investigated the responses of gene regulatory networks containing many genes that have undergone numerical evolution to achieve high fitness due to the adaptive response of only a single target gene; this single target gene responds to changes in external inputs and later returns to basal levels. Despite setting a single target, most genes showed adaptive responses after evolution. Such adaptive dynamics were not due to common motifs within a few genes; even without such motifs, almost all genes showed adaptation, albeit sometimes partial adaptation, in the sense that expression levels did not always return to original levels. The genes split into two groups: genes in the first group exhibited an initial increase in expression and then returned to basal levels, while genes in the second group exhibited the opposite changes in expression. From this model, genes in the first group received positive input from other genes within the first group, but negative input from genes in the second group, and vice versa. Thus, the adaptation dynamics of genes from both groups were consolidated. This cooperative adaptive behavior was commonly observed if the number of genes involved was larger than the order of ten. These results have implications in the collective responses of gene expression networks in microarray measurements of yeast Saccharomyces cerevisiae and the significance to the biological homeostasis of systems with many components., Author Summary Homeostasis is an inherent property of biological systems, which have a general tendency to adapt, i.e., to recover their original state following environmental changes. In cells, this adaptation is mediated by changes in protein expression. Initially, cells respond to environmental changes by altered gene/protein expression; subsequently, the expression of most genes returns to basal levels, albeit not completely, as shown by recent experimental analyses of yeast. Although simple mechanisms for adaptation through network motifs, composed of just a few genes, are well understood, how regulatory networks involving many genes that activate or repress each other can generate adaptive behaviors is unclear. Here, by numerically evolving gene regulatory networks, we obtained a class of genes whose expression dynamics showed adaptation over almost all genes, from which we revealed the general logic underlying such adaptive dynamics with many degrees of freedom, which was not reducible to motifs with a few genes. This adaptation was cooperative, i.e., adaptation of one gene mutually relied upon others' adaptive expressions. Moreover, this collective behavior was robust to noise and mutations. The present study sheds a light on the nature of collective gene expression dynamics allowing for biological homeostasis.

Cellular signal transduction generally involves cascades of post-translational protein modifications that rapidly catalyze changes in protein-DNA interactions and gene expression. High-throughput measurements are improving our ability to study each of these stages individually, but do not capture the connections between them. Here we present an approach for building a network of physical links among these data that can be used to prioritize targets for pharmacological intervention. Our method recovers the critical missing links between proteomic and transcriptional data by relating changes in chromatin accessibility to changes in expression and then uses these links to connect proteomic and transcriptome data. We applied our approach to integrate epigenomic, phosphoproteomic and transcriptome changes induced by the variant III mutation of the epidermal growth factor receptor (EGFRvIII) in a cell line model of glioblastoma multiforme (GBM). To test the relevance of the network, we used small molecules to target highly connected nodes implicated by the network model that were not detected by the experimental data in isolation and we found that a large fraction of these agents alter cell viability. Among these are two compounds, ICG-001, targeting CREB binding protein (CREBBP), and PKF118–310, targeting β-catenin (CTNNB1), which have not been tested previously for effectiveness against GBM. At the level of transcriptional regulation, we used chromatin immunoprecipitation sequencing (ChIP-Seq) to experimentally determine the genome-wide binding locations of p300, a transcriptional co-regulator highly connected in the network. Analysis of p300 target genes suggested its role in tumorigenesis. We propose that this general method, in which experimental measurements are used as constraints for building regulatory networks from the interactome while taking into account noise and missing data, should be applicable to a wide range of high-throughput datasets., National Science Foundation (U.S.) (DB1-0821391), National Institutes of Health (U.S.) (Grant U54-CA112967), National Institutes of Health (U.S.) (Grant R01-GM089903), National Institutes of Health (U.S.) (P30-ES002109)

We computationally determined miRs that are significantly connected to molecular pathways by utilizing gene expression profiles in different cancer types such as glioblastomas, ovarian and breast cancers. Specifically, we assumed that the knowledge of physical interactions between miRs and genes indicated subsets of important miRs (IM) that significantly contributed to the regression of pathway-specific enrichment scores. Despite the different nature of the considered cancer types, we found strongly overlapping sets of IMs. Furthermore, IMs that were important for many pathways were enriched with literature-curated cancer and differentially expressed miRs. Such sets of IMs also coincided well with clusters of miRs that were experimentally indicated in numerous other cancer types. In particular, we focused on an overlapping set of 99 overall important miRs (OIM) that were found in glioblastomas, ovarian and breast cancers simultaneously. Notably, we observed that interactions between OIMs and leading edge genes of differentially expressed pathways were characterized by considerable changes in their expression correlations. Such gains/losses of miR and gene expression correlation indicated miR/gene pairs that may play a causal role in the underlying cancers., Author Summary We assume that a network of physical interactions between miRs and genes allows us to determine miRs that influence the expression of whole pathways in different tumor types. Specifically, we represented each pathway by an enrichment score and an array of miRs counting the number of genes in the pathway a given miR can bind. Despite the different nature of the considered tumor types, we obtained a large set of overlapping miRs using a machine-learning algorithm. Such associated miRs were enriched with literature-curated cancer and differentially expressed miRs and also coincided well with clusters of miRs that were experimentally indicated in numerous other cancer types. Focusing on such sets of miRs we observed that interactions with genes in differentially expressed pathways were characterized by massive gains/losses of expression correlations. Such drastic changes of miR and gene expression correlation indicate miR/gene pairs that may play a causal role in the underlying cancers.

The circadian clock controls many physiological processes in higher plants and causes a large fraction of the genome to be expressed with a 24h rhythm. The transcripts encoding the RNA-binding proteins AtGRP7 (Arabidopsis thaliana Glycine Rich Protein 7) and AtGRP8 oscillate with evening peaks. The circadian clock components CCA1 and LHY negatively affect AtGRP7 expression at the level of transcription. AtGRP7 and AtGRP8, in turn, negatively auto-regulate and reciprocally cross-regulate post-transcriptionally: high protein levels promote the generation of an alternative splice form that is rapidly degraded. This clock-regulated feedback loop has been proposed to act as a molecular slave oscillator in clock output. While mathematical models describing the circadian core oscillator in Arabidopsis thaliana were introduced recently, we propose here the first model of a circadian slave oscillator. We define the slave oscillator in terms of ordinary differential equations and identify the model's parameters by an optimization procedure based on experimental results. The model successfully reproduces the pertinent experimental findings such as waveforms, phases, and half-lives of the time-dependent concentrations. Furthermore, we obtain insights into possible mechanisms underlying the observed experimental dynamics: the negative auto-regulation and reciprocal cross-regulation via alternative splicing could be responsible for the sharply peaking waveforms of the AtGRP7 and AtGRP8 mRNA. Moreover, our results suggest that the AtGRP8 transcript oscillations are subordinated to those of AtGRP7 due to a higher impact of AtGRP7 protein on alternative splicing of its own and of the AtGRP8 pre-mRNA compared to the impact of AtGRP8 protein. Importantly, a bifurcation analysis provides theoretical evidence that the slave oscillator could be a toggle switch, arising from the reciprocal cross-regulation at the post-transcriptional level. In view of this, transcriptional repression of AtGRP7 and AtGRP8 by LHY and CCA1 induces oscillations of the toggle switch, leading to the observed high-amplitude oscillations of AtGRP7 mRNA., Author Summary The circadian clock organizes the day in the life of a plant by causing 24h rhythms in gene expression. For example, the core clockwork of the model plant Arabidopsis thaliana causes the transcripts encoding the RNA-binding proteins AtGRP7 and AtGRP8 to undergo high amplitude oscillations with a peak at the end of the day. AtGRP7 and AtGRP8, in turn, negatively auto-regulate and reciprocally cross-regulate their own expression by causing alternative splicing of their pre-mRNAs, followed by rapid degradation of the alternatively spliced transcripts. This has led to the suggestion that they represent molecular slave oscillators downstream of the core clock. Using a mathematical model we obtain insights into possible mechanisms underlying the experimentally observed dynamics, e.g. a higher impact of AtGRP7 protein compared to the impact of AtGRP8 protein on the alternative splicing explains the experimentally observed phases of their transcript. Previously, components that reciprocally repress their own transcription (double negative loops) have been shown to potentially act as a toggle switch between two states. We provide theoretical evidence that the slave oscillator could be a bistable toggle switch as well, operating at the post-transcriptional level.