Epigenetic factor competition reshapes the EMT landscape.

The emergence of and transitions between distinct phenotypes in isogenic cells can be attributed to the intricate interplay of epigenetic marks, external signals, and gene-regulatory elements. These elements include chromatin remodelers, histone modifiers, transcription factors, and regulatory RNAs. Mathematical models known as gene-regulatory networks (GRNs) are an increasingly important tool to unravel the workings of such complex networks. In such models, epigenetic factors are usually proposed to act on the chromatin regions directly involved in the expression of relevant genes. However, it has been well-established that these factors operate globally and compete with each other for targets genome-wide. Therefore, a perturbation of the activity of a regulator can redistribute epigenetic marks across the genome and modulate the levels of competing regulators. In this paper, we propose a conceptual and mathematical modeling framework that incorporates both local and global competition effects between antagonistic epigenetic regulators, in addition to local transcription factors, and show the counterintuitive consequences of such interactions. We apply our approach to recent experimental findings on the epithelial-mesenchymal transition (EMT). We show that it can explain the puzzling experimental data, as well as provide verifiable predictions.

A novel COVID-19 epidemiological model with explicit susceptible and asymptomatic isolation compartments reveals unexpected consequences of timing social distancing.

Motivated by the current COVID-19 epidemic, this work introduces an epidemiological model in which separate compartments are used for susceptible and asymptomatic "socially distant" populations. Distancing directives are represented by rates of flow into these compartments, as well as by a reduction in contacts that lessens disease transmission. The dynamical behavior of this system is analyzed, under various different rate control strategies, and the sensitivity of the basic reproduction number to various parameters is studied. One of the striking features of this model is the existence of a critical implementation delay (CID) in issuing distancing mandates: while a delay of about two weeks does not have an appreciable effect on the peak number of infections, issuing mandates even slightly after this critical time results in a far greater incidence of infection. Thus, there is a nontrivial but tight "window of opportunity" for commencing social distancing in order to meet the capacity of healthcare resources. However, if one wants to also delay the timing of peak infections - so as to take advantage of potential new therapies and vaccines - action must be taken much faster than the CID. Different relaxation strategies are also simulated, with surprising results. Periodic relaxation policies suggest a schedule which may significantly inhibit peak infective load, but that this schedule is very sensitive to parameter values and the schedule's frequency. Furthermore, we considered the impact of steadily reducing social distancing measures over time. We find that a too-sudden reopening of society may negate the progress achieved under initial distancing guidelines, but the negative effects can be mitigated if the relaxation strategy is carefully designed.

A computational framework for a Lyapunov-enabled analysis of biochemical reaction networks.

Complex molecular biological processes such as transcription and translation, signal transduction, post-translational modification cascades, and metabolic pathways can be described in principle by biochemical reactions that explicitly take into account the sophisticated network of chemical interactions regulating cell life. The ability to deduce the possible qualitative behaviors of such networks from a set of reactions is a central objective and an ongoing challenge in the field of systems biology. Unfortunately, the construction of complete mathematical models is often hindered by a pervasive problem: despite the wealth of qualitative graphical knowledge about network interactions, the form of the governing nonlinearities and/or the values of kinetic constants are hard to uncover experimentally. The kinetics can also change with environmental variations. This work addresses the following question: given a set of reactions and without assuming a particular form for the kinetics, what can we say about the asymptotic behavior of the network? Specifically, it introduces a class of networks that are "structurally (mono) attractive" meaning that they are incapable of exhibiting multiple steady states, oscillation, or chaos by virtue of their reaction graphs. These networks are characterized by the existence of a universal energy-like function called a Robust Lyapunov function (RLF). To find such functions, a finite set of rank-one linear systems is introduced, which form the extremals of a linear convex cone. The problem is then reduced to that of finding a common Lyapunov function for this set of extremals. Based on this characterization, a computational package, Lyapunov-Enabled Analysis of Reaction Networks (LEARN), is provided that constructs such functions or rules out their existence. An extensive study of biochemical networks demonstrates that LEARN offers a new unified framework. Basic motifs, three-body binding, and genetic networks are studied first. The work then focuses on cellular signalling networks including various post-translational modification cascades, phosphotransfer and phosphorelay networks, T-cell kinetic proofreading, and ERK signalling. The Ribosome Flow Model is also studied.

Universal features of epidemic models under social distancing guidelines.

Social distancing as a form of nonpharmaceutical intervention has been enacted in many countries as a form of mitigating the spread of COVID-19. There has been a large interest in mathematical modeling to aid in the prediction of both the total infected population and virus-related deaths, as well as to aid government agencies in decision making. As the virus continues to spread, there are both economic and sociological incentives to minimize time spent with strict distancing mandates enforced, and/or to adopt periodically relaxed distancing protocols, which allow for scheduled economic activity. The main objective of this study is to reduce the disease burden in a population, here measured as the peak of the infected population, while simultaneously minimizing the length of time the population is socially distanced, utilizing both a single period of social distancing as well as periodic relaxation. We derive a linear relationship among the optimal start time and duration of a single interval of social distancing from an approximation of the classic epidemic SIR model. Furthermore, we see a sharp phase transition region in start times for a single pulse of distancing, where the peak of the infected population changes rapidly; notably, this transition occurs well before one would intuitively expect. By numerical investigation of more sophisticated epidemiological models designed specifically to describe the COVID-19 pandemic, we see that all share remarkably similar dynamic characteristics when contact rates are subject to periodic or one-shot changes, and hence lead us to conclude that these features are universal in epidemic models. On the other hand, the nonlinearity of epidemic models leads to non-monotone behavior of the peak of infected population under periodic relaxation of social distancing policies. This observation led us to hypothesize that an additional single interval social distancing at a proper time can significantly decrease the infected peak of periodic policies, and we verified this improvement numerically. While synchronous quarantine and social distancing mandates across populations effectively minimize the spread of an epidemic over the world, relaxation decisions should not be enacted at the same time for different populations.

Oscillatory stimuli differentiate adapting circuit topologies.

Biology emerges from interactions between molecules, which are challenging to elucidate with current techniques. An orthogonal approach is to probe for 'response signatures' that identify specific circuit motifs. For example, bistability, hysteresis, or irreversibility are used to detect positive feedback loops. For adapting systems, such signatures are not known. Only two circuit motifs generate adaptation: negative feedback loops (NFLs) and incoherent feed-forward loops (IFFLs). On the basis of computational testing and mathematical proofs, we propose differential signatures: in response to oscillatory stimulation, NFLs but not IFFLs show refractory-period stabilization (robustness to changes in stimulus duration) or period skipping. Applying this approach to yeast, we identified the circuit dominating cell cycle timing. In Caenorhabditis elegans AWA neurons, which are crucial for chemotaxis, we uncovered a Ca2+ NFL leading to adaptation that would be difficult to find by other means. These response signatures allow direct access to the outlines of the wiring diagrams of adapting systems.

Integrating transcriptomics and bulk time course data into a mathematical framework to describe and predict therapeutic resistance in cancer.

A significant challenge in the field of biomedicine is the development of methods to integrate the multitude of dispersed data sets into comprehensive frameworks to be used to generate optimal clinical decisions. Recent technological advances in single cell analysis allow for high-dimensional molecular characterization of cells and populations, but to date, few mathematical models have attempted to integrate measurements from the single cell scale with other types of longitudinal data. Here, we present a framework that actionizes static outputs from a machine learning model and leverages these as measurements of state variables in a dynamic model of treatment response. We apply this framework to breast cancer cells to integrate single cell transcriptomic data with longitudinal bulk cell population (bulk time course) data. We demonstrate that the explicit inclusion of the phenotypic composition estimate, derived from single cell RNA-sequencing data (scRNA-seq), improves accuracy in the prediction of new treatments with a concordance correlation coefficient (CCC) of 0.92 compared to a prediction accuracy of CCC = 0.64 when fitting on longitudinal bulk cell population data alone. To our knowledge, this is the first work that explicitly integrates single cell clonally-resolved transcriptome datasets with bulk time-course data to jointly calibrate a mathematical model of drug resistance dynamics. We anticipate this approach to be a first step that demonstrates the feasibility of incorporating multiple data types into mathematical models to develop optimized treatment regimens from data.

Inferring reaction network structure from single-cell, multiplex data, using toric systems theory.

The goal of many single-cell studies on eukaryotic cells is to gain insight into the biochemical reactions that control cell fate and state. In this paper we introduce the concept of Effective Stoichiometric Spaces (ESS) to guide the reconstruction of biochemical networks from multiplexed, fixed time-point, single-cell data. In contrast to methods based solely on statistical models of data, the ESS method leverages the power of the geometric theory of toric varieties to begin unraveling the structure of chemical reaction networks (CRN). This application of toric theory enables a data-driven mapping of covariance relationships in single-cell measurements into stoichiometric information, one in which each cell subpopulation has its associated ESS interpreted in terms of CRN theory. In the development of ESS we reframe certain aspects of the theory of CRN to better match data analysis. As an application of our approach we process cytomery- and image-based single-cell datasets and identify differences in cells treated with kinase inhibitors. Our approach is directly applicable to data acquired using readily accessible experimental methods such as Fluorescence Activated Cell Sorting (FACS) and multiplex immunofluorescence.

Multi-modality in gene regulatory networks with slow promoter kinetics.

Phenotypical variability in the absence of genetic variation often reflects complex energetic landscapes associated with underlying gene regulatory networks (GRNs). In this view, different phenotypes are associated with alternative states of complex nonlinear systems: stable attractors in deterministic models or modes of stationary distributions in stochastic descriptions. We provide theoretical and practical characterizations of these landscapes, specifically focusing on stochastic Slow Promoter Kinetics (SPK), a time scale relevant when transcription factor binding and unbinding are affected by epigenetic processes like DNA methylation and chromatin remodeling. In this case, largely unexplored except for numerical simulations, adiabatic approximations of promoter kinetics are not appropriate. In contrast to the existing literature, we provide rigorous analytic characterizations of multiple modes. A general formal approach gives insight into the influence of parameters and the prediction of how changes in GRN wiring, for example through mutations or artificial interventions, impact the possible number, location, and likelihood of alternative states. We adapt tools from the mathematical field of singular perturbation theory to represent stationary distributions of Chemical Master Equations for GRNs as mixtures of Poisson distributions and obtain explicit formulas for the locations and probabilities of metastable states as a function of the parameters describing the system. As illustrations, the theory is used to tease out the role of cooperative binding in stochastic models in comparison to deterministic models, and applications are given to various model systems, such as toggle switches in isolation or in communicating populations, a synthetic oscillator, and a trans-differentiation network.

Evaluating optimal therapy robustness by virtual expansion of a sample population, with a case study in cancer immunotherapy.

Cancer is a highly heterogeneous disease, exhibiting spatial and temporal variations that pose challenges for designing robust therapies. Here, we propose the VEPART (Virtual Expansion of Populations for Analyzing Robustness of Therapies) technique as a platform that integrates experimental data, mathematical modeling, and statistical analyses for identifying robust optimal treatment protocols. VEPART begins with time course experimental data for a sample population, and a mathematical model fit to aggregate data from that sample population. Using nonparametric statistics, the sample population is amplified and used to create a large number of virtual populations. At the final step of VEPART, robustness is assessed by identifying and analyzing the optimal therapy (perhaps restricted to a set of clinically realizable protocols) across each virtual population. As proof of concept, we have applied the VEPART method to study the robustness of treatment response in a mouse model of melanoma subject to treatment with immunostimulatory oncolytic viruses and dendritic cell vaccines. Our analysis (i) showed that every scheduling variant of the experimentally used treatment protocol is fragile (nonrobust) and (ii) discovered an alternative region of dosing space (lower oncolytic virus dose, higher dendritic cell dose) for which a robust optimal protocol exists.

Subharmonics and Chaos in Simple Periodically Forced Biomolecular Models.

This article uncovers a remarkable behavior in two biochemical systems that commonly appear as components of signal transduction pathways in systems biology. These systems have globally attracting steady states when unforced, so they might have been considered uninteresting from a dynamical standpoint. However, when subject to a periodic excitation, strange attractors arise via a period-doubling cascade. Quantitative analyses of the corresponding discrete chaotic trajectories are conducted numerically by computing largest Lyapunov exponents, power spectra, and autocorrelation functions. To gain insight into the geometry of the strange attractors, the phase portraits of the corresponding iterated maps are interpreted as scatter plots for which marginal distributions are additionally evaluated. The lack of entrainment to external oscillations, in even the simplest biochemical networks, represents a level of additional complexity in molecular biology, which has previously been insufficiently recognized but is plausibly biologically important.

Dynamic compensation, parameter identifiability, and equivariances.

A recent paper by Karin et al. introduced a mathematical notion called dynamical compensation (DC) of biological circuits. DC was shown to play an important role in glucose homeostasis as well as other key physiological regulatory mechanisms. Karin et al. went on to provide a sufficient condition to test whether a given system has the DC property. Here, we show how DC can be formulated in terms of a well-known concept in systems biology, statistics, and control theory-that of parameter structural non-identifiability. Viewing DC as a parameter identification problem enables one to take advantage of powerful theoretical and computational tools to test a system for DC. We obtain as a special case the sufficient criterion discussed by Karin et al. We also draw connections to system equivalence and to the fold-change detection property.

Reduction of multiscale stochastic biochemical reaction networks using exact moment derivation.

Biochemical reaction networks (BRNs) in a cell frequently consist of reactions with disparate timescales. The stochastic simulations of such multiscale BRNs are prohibitively slow due to high computational cost for the simulations of fast reactions. One way to resolve this problem uses the fact that fast species regulated by fast reactions quickly equilibrate to their stationary distribution while slow species are unlikely to be changed. Thus, on a slow timescale, fast species can be replaced by their quasi-steady state (QSS): their stationary conditional expectation values for given slow species. As the QSS are determined solely by the state of slow species, such replacement leads to a reduced model, where fast species are eliminated. However, it is challenging to derive the QSS in the presence of nonlinear reactions. While various approximation schemes for the QSS have been developed, they often lead to considerable errors. Here, we propose two classes of multiscale BRNs which can be reduced by deriving an exact QSS rather than approximations. Specifically, if fast species constitute either a feedforward network or a complex balanced network, the reduced model based on the exact QSS can be derived. Such BRNs are frequently observed in a cell as the feedforward network is one of fundamental motifs of gene or protein regulatory networks. Furthermore, complex balanced networks also include various types of fast reversible bindings such as bindings between transcriptional factors and gene regulatory sites. The reduced models based on exact QSS, which can be calculated by the computational packages provided in this work, accurately approximate the slow scale dynamics of the original full model with much lower computational cost.

Computation-Guided Design of a Stimulus-Responsive Multienzyme Supramolecular Assembly.

The construction of stimulus-responsive supramolecular complexes of metabolic pathway enzymes, inspired by natural multienzyme assemblies (metabolons), provides an attractive avenue for efficient and spatiotemporally controllable one-pot biotransformations. We have constructed a phosphorylation- and optically responsive metabolon for the biodegradation of the environmental pollutant 1,2,3-trichloropropane.

Quorum-Sensing Synchronization of Synthetic Toggle Switches: A Design Based on Monotone Dynamical Systems Theory.

Synthetic constructs in biotechnology, biocomputing, and modern gene therapy interventions are often based on plasmids or transfected circuits which implement some form of "on-off" switch. For example, the expression of a protein used for therapeutic purposes might be triggered by the recognition of a specific combination of inducers (e.g., antigens), and memory of this event should be maintained across a cell population until a specific stimulus commands a coordinated shut-off. The robustness of such a design is hampered by molecular ("intrinsic") or environmental ("extrinsic") noise, which may lead to spontaneous changes of state in a subset of the population and is reflected in the bimodality of protein expression, as measured for example using flow cytometry. In this context, a "majority-vote" correction circuit, which brings deviant cells back into the required state, is highly desirable, and quorum-sensing has been suggested as a way for cells to broadcast their states to the population as a whole so as to facilitate consensus. In this paper, we propose what we believe is the first such a design that has mathematically guaranteed properties of stability and auto-correction under certain conditions. Our approach is guided by concepts and theory from the field of "monotone" dynamical systems developed by M. Hirsch, H. Smith, and others. We benchmark our design by comparing it to an existing design which has been the subject of experimental and theoretical studies, illustrating its superiority in stability and self-correction of synchronization errors. Our stability analysis, based on dynamical systems theory, guarantees global convergence to steady states, ruling out unpredictable ("chaotic") behaviors and even sustained oscillations in the limit of convergence. These results are valid no matter what are the values of parameters, and are based only on the wiring diagram. The theory is complemented by extensive computational bifurcation analysis, performed for a biochemically-detailed and biologically-relevant model that we developed. Another novel feature of our approach is that our theorems on exponential stability of steady states for homogeneous or mixed populations are valid independently of the number N of cells in the population, which is usually very large (N ≫ 1) and unknown. We prove that the exponential stability depends on relative proportions of each type of state only. While monotone systems theory has been used previously for systems biology analysis, the current work illustrates its power for synthetic biology design, and thus has wider significance well beyond the application to the important problem of coordination of toggle switches.

Minimization of thermodynamic costs in cancer cell invasion.

Metastasis, the truly lethal aspect of cancer, occurs when metastatic cancer cells in a tumor break through the basement membrane and penetrate the extracellular matrix. We show that MDA-MB-231 metastatic breast cancer cells cooperatively invade a 3D collagen matrix while following a glucose gradient. The invasion front of the cells is a dynamic one, with different cells assuming the lead on a time scale of 70 h. The front cell leadership is dynamic presumably because of metabolic costs associated with a long-range strain field that precedes the invading cell front, which we have imaged using confocal imaging and marker beads imbedded in the collagen matrix. We suggest this could be a quantitative assay for an invasive phenotype tracking a glucose gradient and show that the invading cells act in a cooperative manner by exchanging leaders in the invading front.

A 'resource allocator' for transcription based on a highly fragmented T7 RNA polymerase.

Synthetic genetic systems share resources with the host, including machinery for transcription and translation. Phage RNA polymerases (RNAPs) decouple transcription from the host and generate high expression. However, they can exhibit toxicity and lack accessory proteins (σ factors and activators) that enable switching between different promoters and modulation of activity. Here, we show that T7 RNAP (883 amino acids) can be divided into four fragments that have to be co-expressed to function. The DNA-binding loop is encoded in a C-terminal 285-aa 'σ fragment', and fragments with different specificity can direct the remaining 601-aa 'core fragment' to different promoters. Using these parts, we have built a resource allocator that sets the core fragment concentration, which is then shared by multiple σ fragments. Adjusting the concentration of the core fragment sets the maximum transcriptional capacity available to a synthetic system. Further, positive and negative regulation is implemented using a 67-aa N-terminal 'α fragment' and a null (inactivated) σ fragment, respectively. The α fragment can be fused to recombinant proteins to make promoters responsive to their levels. These parts provide a toolbox to allocate transcriptional resources via different schemes, which we demonstrate by building a system which adjusts promoter activity to compensate for the difference in copy number of two plasmids.

Mechanism-independent method for predicting response to multidrug combinations in bacteria.

Drugs are commonly used in combinations larger than two for treating bacterial infection. However, it is generally impossible to infer directly from the effects of individual drugs the net effect of a multidrug combination. Here we develop a mechanism-independent method for predicting the microbial growth response to combinations of more than two drugs. Performing experiments in both Gram-negative (Escherichia coli) and Gram-positive (Staphylococcus aureus) bacteria, we demonstrate that for a wide range of drugs, the bacterial responses to drug pairs are sufficient to infer the effects of larger drug combinations. To experimentally establish the broad applicability of the method, we use drug combinations comprising protein synthesis inhibitors (macrolides, aminoglycosides, tetracyclines, lincosamides, and chloramphenicol), DNA synthesis inhibitors (fluoroquinolones and quinolones), folic acid synthesis inhibitors (sulfonamides and diaminopyrimidines), cell wall synthesis inhibitors, polypeptide antibiotics, preservatives, and analgesics. Moreover, we show that the microbial responses to these drug combinations can be predicted using a simple formula that should be widely applicable in pharmacology. These findings offer a powerful, readily accessible method for the rational design of candidate therapies using combinations of more than two drugs. In addition, the accurate predictions of this framework raise the question of whether the multidrug response in bacteria obeys statistical, rather than chemical, laws for combinations larger than two.

The energy costs of insulators in biochemical networks.

Complex networks of biochemical reactions, such as intracellular protein signaling pathways and genetic networks, are often conceptualized in terms of modules--semiindependent collections of components that perform a well-defined function and which may be incorporated in multiple pathways. However, due to sequestration of molecular messengers during interactions and other effects, collectively referred to as retroactivity, real biochemical systems do not exhibit perfect modularity. Biochemical signaling pathways can be insulated from impedance and competition effects, which inhibit modularity, through enzymatic futile cycles that consume energy, typically in the form of ATP. We hypothesize that better insulation necessarily requires higher energy consumption. We test this hypothesis through a combined theoretical and computational analysis of a simplified physical model of covalent cycles, using two innovative measures of insulation, as well as a possible new way to characterize optimal insulation through the balancing of these two measures in a Pareto sense. Our results indicate that indeed better insulation requires more energy. While insulation may facilitate evolution by enabling a modular plug-and-play interconnection architecture, allowing for the creation of new behaviors by adding targets to existing pathways, our work suggests that this potential benefit must be balanced against the metabolic costs of insulation necessarily incurred in not affecting the behavior of existing processes.

A characterization of scale invariant responses in enzymatic networks.

An ubiquitous property of biological sensory systems is adaptation: a step increase in stimulus triggers an initial change in a biochemical or physiological response, followed by a more gradual relaxation toward a basal, pre-stimulus level. Adaptation helps maintain essential variables within acceptable bounds and allows organisms to readjust themselves to an optimum and non-saturating sensitivity range when faced with a prolonged change in their environment. Recently, it was shown theoretically and experimentally that many adapting systems, both at the organism and single-cell level, enjoy a remarkable additional feature: scale invariance, meaning that the initial, transient behavior remains (approximately) the same even when the background signal level is scaled. In this work, we set out to investigate under what conditions a broadly used model of biochemical enzymatic networks will exhibit scale-invariant behavior. An exhaustive computational study led us to discover a new property of surprising simplicity and generality, uniform linearizations with fast output (ULFO), whose validity we show is both necessary and sufficient for scale invariance of three-node enzymatic networks (and sufficient for any number of nodes). Based on this study, we go on to develop a mathematical explanation of how ULFO results in scale invariance. Our work provides a surprisingly consistent, simple, and general framework for understanding this phenomenon, and results in concrete experimental predictions.

Logarithmic Lipschitz norms and diffusion-induced instability.

This paper proves that contractive ordinary differential equation systems remain contractive when diffusion is added. Thus, diffusive instabilities, in the sense of the Turing phenomenon, cannot arise for such systems. An important biochemical system is shown to satisfy the required conditions.