Quantifying robustness of the gap gene network.

Early development of Drosophila melanogaster (fruit fly) facilitated by the gap gene network has been shown to be incredibly robust, and the same patterns emerge even when the process is seriously disrupted. We investigate this robustness using a previously developed computational framework called DSGRN (Dynamic Signatures Generated by Regulatory Networks). Our mathematical innovations include the conceptual extension of this established modeling technique to enable modeling of spatially monotone environmental effects, as well as the development of a collection of graph theoretic robustness scores for network models. This allows us to rank order the robustness of network models of cellular systems where each cell contains the same genetic network topology but operates under a parameter regime that changes continuously from cell to cell. We demonstrate the power of this method by comparing the robustness of two previously introduced network models of gap gene expression along the anterior-posterior axis of the fruit fly embryo, both to each other and to a random sample of networks with same number of nodes and edges. We observe that there is a substantial difference in robustness scores between the two models. Our biological insight is that random network topologies are in general capable of reproducing complex patterns of expression, but that using measures of robustness to rank order networks permits a large reduction in hypothesis space for highly conserved systems such as developmental networks.

Experimental guidance for discovering genetic networks through hypothesis reduction on time series.

Large programs of dynamic gene expression, like cell cyles and circadian rhythms, are controlled by a relatively small "core" network of transcription factors and post-translational modifiers, working in concerted mutual regulation. Recent work suggests that system-independent, quantitative features of the dynamics of gene expression can be used to identify core regulators. We introduce an approach of iterative network hypothesis reduction from time-series data in which increasingly complex features of the dynamic expression of individual, pairs, and entire collections of genes are used to infer functional network models that can produce the observed transcriptional program. The culmination of our work is a computational pipeline, Iterative Network Hypothesis Reduction from Temporal Dynamics (Inherent dynamics pipeline), that provides a priority listing of targets for genetic perturbation to experimentally infer network structure. We demonstrate the capability of this integrated computational pipeline on synthetic and yeast cell-cycle data.

Ribosome Abundance Control in Prokaryotes.

Cell growth is an essential phenotype of any unicellular organism and it crucially depends on precise control of protein synthesis. We construct a model of the feedback mechanisms that regulate abundance of ribosomes in E. coli, a prototypical prokaryotic organism. Since ribosomes are needed to produce more ribosomes, the model includes a positive feedback loop central to the control of cell growth. Our analysis of the model shows that there can be only two coexisting equilibrium states across all 23 parameters. This precludes the existence of hysteresis, suggesting that the ribosome abundance changes continuously with parameters. These states are related by a transcritical bifurcation, and we provide an analytic formula for parameters that admit either state.

A Nonlinear Delay Model for Metabolic Oscillations in Yeast Cells.

We introduce two time-delay models of metabolic oscillations in yeast cells. Our model tests a hypothesis that the oscillations occur as multiple pathways share a limited resource which we equate to the number of available ribosomes. We initially explore a single-protein model with a constraint equation governing the total resource available to the cell. The model is then extended to include three proteins that share a resource pool. Three approaches are considered at constant delay to numerically detect oscillations. First, we use a spectral element method to approximate the system as a discrete map and evaluate the stability of the linearized system about its equilibria by examining its eigenvalues. For the second method, we plot amplitudes of the simulation trajectories in 2D projections of the parameter space. We use a history function that is consistent with published experimental results to obtain metabolic oscillations. Finally, the spectral element method is used to convert the system to a boundary value problem whose solutions correspond to approximate periodic solutions of the system. Our results show that certain combinations of total resource available and the time delay, lead to oscillations. We observe that an oscillation region in the parameter space is between regions admitting steady states that correspond to zero and constant production. Similar behavior is found with the three-protein model where all proteins require the same production time. However, a shift in the protein production rates peaks occurs for low available resource suggesting that our model captures the shared resource pool dynamics.

Genetic networks encode secrets of their past.

Research shows that gene duplication followed by either repurposing or removal of duplicated genes is an important contributor to evolution of gene and protein interaction networks. We aim to identify which characteristics of a network can arise through this process, and which must have been produced in a different way. To model the network evolution, we postulate vertex duplication and edge deletion as evolutionary operations on graphs. Using the novel concept of an ancestrally distinguished subgraph, we show how features of present-day networks require certain features of their ancestors. In particular, ancestrally distinguished subgraphs cannot be introduced by vertex duplication. Additionally, if vertex duplication and edge deletion are the only evolutionary mechanisms, then a graph's ancestrally distinguished subgraphs must be contained in all of the graph's ancestors. We analyze two experimentally derived genetic networks and show that our results accurately predict lack of large ancestrally distinguished subgraphs, despite this feature being statistically improbable in associated random networks. This observation is consistent with the hypothesis that these networks evolved primarily via vertex duplication. The tools we provide open the door for analyzing ancestral networks using current networks. Our results apply to edge-labeled (e.g. signed) graphs which are either undirected or directed.

Rational design of complex phenotype via network models.

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.

Positive feedback drives a secondary nonlinear product burst during a biphasic DNA amplification reaction.

Isothermal DNA amplification reactions are used in a broad variety of applications, from diagnostic assays to DNA circuits, with greater speed and less complexity than established PCR technologies. We recently reported a unique, high gain, biphasic isothermal DNA amplification reaction, called the Ultrasensitive DNA Amplification Reaction (UDAR). Here we present a detailed analysis of the UDAR reaction pathways that initiates with a first phase followed by a nonlinear product burst, which is caused by an autocatalytic secondary reaction. The experimental reaction output was reproduced using an ordinary differential equation model based on detailed reaction mechanisms. This model provides insight on the relative importance of each reaction mechanism during both phases, which could aid in the design of product output during DNA amplification reactions.

Modeling Transport Regulation in Gene Regulatory Networks.

A gene regulatory network summarizes the interactions between a set of genes and regulatory gene products. These interactions include transcriptional regulation, protein activity regulation, and regulation of the transport of proteins between cellular compartments. DSGRN is a network modeling approach that builds on traditions of discrete-time Boolean models and continuous-time switching system models. When all interactions are transcriptional, DSGRN uses a combinatorial approximation to describe the entire range of dynamics that is compatible with network structure. Here we present an extension of the DGSRN approach to transport regulation across a boundary between compartments, such as a cellular membrane. We illustrate our approach by searching a model of the p53-Mdm2 network for the potential to admit two experimentally observed distinct stable periodic cycles.

Operon dynamics with state dependent transcription and/or translation delays.

Transcription and translation retrieve and operationalize gene encoded information in cells. These processes are not instantaneous and incur significant delays. In this paper we study Goodwin models of both inducible and repressible operons with state-dependent delays. The paper provides justification and derivation of the model, detailed analysis of the appropriate setting of the corresponding dynamical system, and extensive numerical analysis of its dynamics. Comparison with constant delay models shows significant differences in dynamics that include existence of stable periodic orbits in inducible systems and multistability in repressible systems. A combination of parameter space exploration, numerics, analysis of steady state linearization and bifurcation theory indicates the likely presence of Shilnikov-type homoclinic bifurcations in the repressible operon model.

Multistability in the epithelial-mesenchymal transition network.

BACKGROUND: The transitions between epithelial (E) and mesenchymal (M) cell phenotypes are essential in many biological processes like tissue development and cancer metastasis. Previous studies, both modeling and experimental, suggested that in addition to E and M states, the network responsible for these phenotypes exhibits intermediate phenotypes between E and M states. The number and importance of such states is subject to intense discussion in the epithelial-mesenchymal transition (EMT) community. RESULTS: Previous modeling efforts used traditional bifurcation analysis to explore the number of the steady states that correspond to E, M and intermediate states by varying one or two parameters at a time. Since the system has dozens of parameters that are largely unknown, it remains a challenging problem to fully describe the potential set of states and their relationship across all parameters. We use the computational tool DSGRN (Dynamic Signatures Generated by Regulatory Networks) to explore the intermediate states of an EMT model network by computing summaries of the dynamics across all of parameter space. We find that the only attractors in the system are equilibria, that E and M states dominate across parameter space, but that bistability and multistability are common. Even at extreme levels of some of the known inducers of the transition, there is a certain proportion of the parameter space at which an E or an M state co-exists with other stable steady states. CONCLUSIONS: Our results suggest that the multistability is broadly present in the EMT network across parameters and thus response of cells to signals may strongly depend on the particular cell line and genetic background.

Using extremal events to characterize noisy time series.

Experimental time series provide an informative window into the underlying dynamical system, and the timing of the extrema of a time series (or its derivative) contains information about its structure. However, the time series often contain significant measurement errors. We describe a method for characterizing a time series for any assumed level of measurement error [Formula: see text] by a sequence of intervals, each of which is guaranteed to contain an extremum for any function that [Formula: see text]-approximates the time series. Based on the merge tree of a continuous function, we define a new object called the normalized branch decomposition, which allows us to compute intervals for any level [Formula: see text]. We show that there is a well-defined total order on these intervals for a single time series, and that it is naturally extended to a partial order across a collection of time series comprising a dataset. We use the order of the extracted intervals in two applications. First, the partial order describing a single dataset can be used to pattern match against switching model output (Cummins et al. in SIAM J Appl Dyn Syst 17(2):1589-1616, 2018), which allows the rejection of a network model. Second, the comparison between graph distances of the partial orders of different datasets can be used to quantify similarity between biological replicates.

Role of resource allocation and transport in emergence of cross-feeding in microbial consortia.

Microbial communities that implement mutual cross-feeding are commonly observed in nature and with synthetic constructs in laboratory experiments. A mathematical model of competition in a chemostat is developed to investigate the role that resource allocation and transport of metabolites play in cooperation. The model contains four cell types that differ by whether they produce two, one, or none of two essential metabolites. Producing cell types may export these resources into the environment, and those that do not produce both metabolites must import the missing resource. The contribution to the emergence of a collaborative consortium of single resource producers from the transport rate of these metabolites and the type of transport used by the cell (active vs. passive) is studied. Multiple instances of bi-stability and tri-stability are observed, and the effect of the initial concentration of a non-cooperative cheater cell type on the final outcome of the competition is examined. When the cost of producing metabolites is introduced into the model, significant changes to the outcome of the competition are observed, including coexistence of multiple cell types.

Identifying robust hysteresis in networks.

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.

Comparison of Combinatorial Signatures of Global Network Dynamics Generated by Two Classes of ODE Models.

Modeling the dynamics of biological networks introduces many challenges, among them the lack of first principle models, the size of the networks, and difficulties with parameterization. Discrete time Boolean networks and related continuous time switching systems provide a computationally accessible way to translate the structure of the network to predictions about the dynamics. Recent work has shown that the parameterized dynamics of switching systems can be captured by a combinatorial object, called a Dynamic Signatures Generated by Regulatory Networks (DSGRN) database, that consists of a parameter graph characterizing a finite parameter space decomposition, whose nodes are assigned a Morse graph that captures global dynamics for all corresponding parameters. We show that for a given network there is a way to associate the same type of object by considering a continuous time ODE system with a continuous right-hand side, which we call an L-system. The main goal of this paper is to compare the two DSGRN databases for the same network. Since the L-systems can be thought of as perturbations (not necessarily small) of the switching systems, our results address the correspondence between global parameterized dynamics of switching systems and their perturbations. We show that, at corresponding parameters, there is an order preserving map from the Morse graph of the switching system to that of the L-system that is surjective on the set of attractors and bijective on the set of fixed-point attractors. We provide important examples showing why this correspondence cannot be strengthened.

Competitive resource allocation to metabolic pathways contributes to overflow metabolisms and emergent properties in cross-feeding microbial consortia.

Resource scarcity is a common stress in nature and has a major impact on microbial physiology. This review highlights microbial acclimations to resource scarcity, focusing on resource investment strategies for chemoheterotrophs from the molecular level to the pathway level. Competitive resource allocation strategies often lead to a phenotype known as overflow metabolism; the resulting overflow byproducts can stabilize cooperative interactions in microbial communities and can lead to cross-feeding consortia. These consortia can exhibit emergent properties such as enhanced resource usage and biomass productivity. The literature distilled here draws parallels between in silico and laboratory studies and ties them together with ecological theories to better understand microbial stress responses and mutualistic consortia functioning.

Model Rejection and Parameter Reduction via Time Series.

We show how a graph algorithm for finding matching labeled paths in pairs of labeled directed graphs can be used to perform model invalidation for a class of dynamical systems including regulatory network models of relevance to systems biology. In particular, given a partial order of events describing local minima and local maxima of observed quantities from experimental time series data, we produce a labeled directed graph we call the pattern graph for which every path from root to leaf corresponds to a plausible sequence of events. We then consider the regulatory network model, which can itself be rendered into a labeled directed graph we call the search graph via techniques previously developed in computational dynamics. Labels on the pattern graph correspond to experimentally observed events, while labels on the search graph correspond to mathematical facts about the model. We give a theoretical guarantee that failing to find a match invalidates the model. As an application we consider gene regulatory models for the yeast S. cerevisiae.

Global dynamics for switching systems and their extensions by linear differential equations.

Switching systems use piecewise constant nonlinearities to model gene regulatory networks. This choice provides advantages in the analysis of behavior and allows the global description of dynamics in terms of Morse graphs associated to nodes of a parameter graph. The parameter graph captures spatial characteristics of a decomposition of parameter space into domains with identical Morse graphs. However, there are many cellular processes that do not exhibit threshold-like behavior and thus are not well described by a switching system. We consider a class of extensions of switching systems formed by a mixture of switching interactions and chains of variables governed by linear differential equations. We show that the parameter graphs associated to the switching system and any of its extensions are identical. For each parameter graph node, there is an order-preserving map from the Morse graph of the switching system to the Morse graph of any of its extensions. We provide counterexamples that show why possible stronger relationships between the Morse graphs are not valid.

A Mechanistic Model for Cooperative Behavior of Co-transcribing RNA Polymerases.

In fast-transcribing prokaryotic genes, such as an rrn gene in Escherichia coli, many RNA polymerases (RNAPs) transcribe the DNA simultaneously. Active elongation of RNAPs is often interrupted by pauses, which has been observed to cause RNAP traffic jams; yet some studies indicate that elongation seems to be faster in the presence of multiple RNAPs than elongation by a single RNAP. We propose that an interaction between RNAPs via the torque produced by RNAP motion on helically twisted DNA can explain this apparent paradox. We have incorporated the torque mechanism into a stochastic model and simulated transcription both with and without torque. Simulation results illustrate that the torque causes shorter pause durations and fewer collisions between polymerases. Our results suggest that the torsional interaction of RNAPs is an important mechanism in maintaining fast transcription times, and that transcription should be viewed as a cooperative group effort by multiple polymerases.

Convergence Properties of Posttranslationally Modified Protein-Protein Switching Networks with Fast Decay Rates.

A significant conceptual difficulty in the use of switching systems to model regulatory networks is the presence of so-called "black walls," co-dimension 1 regions of phase space with a vector field pointing inward on both sides of the hyperplane. Black walls result from the existence of direct negative self-regulation in the system. One biologically inspired way of removing black walls is the introduction of intermediate variables that mediate the negative self-regulation. In this paper, we study such a perturbation. We replace a switching system with a higher-dimensional switching system with rapidly decaying intermediate proteins, and compare the dynamics between the two systems. We find that the while the individual solutions of the original system can be approximated for a finite time by solutions of a sufficiently close perturbed system, there are always solutions that are not well approximated for any fixed perturbation. We also study a particular example, where global basins of attraction of the perturbed system have a strikingly different form than those of the original system. We perform this analysis using techniques that are adapted to dealing with non-smooth systems.

Cyclic Feedback Systems with Quorum Sensing Coupling.

Synchronization and desynchronization is of great interest in the study of circadian rhythms, metabolic oscillations and time-dependent cell aggregate behaviors. Several recent studies examine synchronization and other dynamics in models of repressilators coupled by a quorum sensing mechanism that uses a diffusive signal. Their numerical simulations have shown the complexity of the collective behavior depends sensitively on which protein upregulates diffusive signal. In this paper, we rigorously prove that the collective dynamics indeed strongly depends on how the signaling network integrates into the repressilator network. In fact we prove a general result for a class of negative cyclic feedback systems with signaling of which the repressilator is but one example. We show that if the feedback along the signaling loop is also negative, the resulting negative feedback, negative signaling (Nf-Ns) system admits either unique stable equilibrium, or a stable oscillation. When a positive signaling feedback is included, the system is no longer (Nf-Ns) and numerically exhibits multistable dynamics (Ullner et al. in Phys Rev Lett 99:148103, 2007; Phys Rev E 78:031904, 2008). We demonstrate that this multistability emerges through saddle node bifurcations of a sole cubic curve-as in generic bistable models.