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.

Solving Immunology?

Emergent responses of the immune system result from the integration of molecular and cellular networks over time and across multiple organs. High-content and high-throughput analysis technologies, concomitantly with data-driven and mechanistic modeling, hold promise for the systematic interrogation of these complex pathways. However, connecting genetic variation and molecular mechanisms to individual phenotypes and health outcomes has proven elusive. Gaps remain in data, and disagreements persist about the value of mechanistic modeling for immunology. Here, we present the perspectives that emerged from the National Institute of Allergy and Infectious Disease (NIAID) workshop 'Complex Systems Science, Modeling and Immunity' and subsequent discussions regarding the potential synergy of high-throughput data acquisition, data-driven modeling, and mechanistic modeling to define new mechanisms of immunological disease and to accelerate the translation of these insights into therapies.

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.

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.

Transcriptional control of human p53-regulated genes.

The p53 protein regulates the transcription of many different genes in response to a wide variety of stress signals. Following DNA damage, p53 regulates key processes, including DNA repair, cell-cycle arrest, senescence and apoptosis, in order to suppress cancer. This Analysis article provides an overview of the current knowledge of p53-regulated genes in these pathways and others, and the mechanisms of their regulation. In addition, we present the most comprehensive list so far of human p53-regulated genes and their experimentally validated, functional binding sites that confer p53 regulation.

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.

Discriminating direct and indirect connectivities in biological networks.

Reverse engineering of biological pathways involves an iterative process between experiments, data processing, and theoretical analysis. Despite concurrent advances in quality and quantity of data as well as computing resources and algorithms, difficulties in deciphering direct and indirect network connections are prevalent. Here, we adopt the notions of abstraction, emulation, benchmarking, and validation in the context of discovering features specific to this family of connectivities. After subjecting benchmark synthetic circuits to perturbations, we inferred the network connections using a combination of nonparametric single-cell data resampling and modular response analysis. Intriguingly, we discovered that recovered weights of specific network edges undergo divergent shifts under differential perturbations, and that the particular behavior is markedly different between topologies. Our results point to a conceptual advance for reverse engineering beyond weight inference. Investigating topological changes under differential perturbations may address the longstanding problem of discriminating direct and indirect connectivities in biological networks.

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.

Non-monotonic Response to Monotonic Stimulus: Regulation of Glyoxylate Shunt Gene-Expression Dynamics in Mycobacterium tuberculosis.

Understanding how dynamical responses of biological networks are constrained by underlying network topology is one of the fundamental goals of systems biology. Here we employ monotone systems theory to formulate a theorem stating necessary conditions for non-monotonic time-response of a biochemical network to a monotonic stimulus. We apply this theorem to analyze the non-monotonic dynamics of the σB-regulated glyoxylate shunt gene expression in Mycobacterium tuberculosis cells exposed to hypoxia. We first demonstrate that the known network structure is inconsistent with observed dynamics. To resolve this inconsistency we employ the formulated theorem, modeling simulations and optimization along with follow-up dynamic experimental measurements. We show a requirement for post-translational modulation of σB activity in order to reconcile the network dynamics with its topology. The results of this analysis make testable experimental predictions and demonstrate wider applicability of the developed methodology to a wide class of biological systems.

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.