Simplifying biochemical models with intermediate species

Connectivity of Parameter Regions of Multistationarity for Multisite Phosphorylation Networks

The parameter region of multistationarity of a reaction network contains all the parameters for which the associated dynamical system exhibits multiple steady states. Describing this region is challenging and remains an active area of research. In this paper, we concentrate on two biologically relevant families of reaction networks that model multisite phosphorylation and dephosphorylation of a substrate at n sites. For small values of n, it had previously been shown that the parameter region of multistationarity is connected. Here, we extend these results and provide a proof that applies to all values of n. Our techniques are based on the study of the critical polynomial associated with these reaction networks together with polyhedral geometric conditions of the signed support of this polynomial.

Multistationarity questions in reduced versus extended biochemical networks

We address several questions in reduced versus extended networks via the elimination or addition of intermediate complexes in the framework of chemical reaction networks with mass-action kinetics. We clarify and extend advances in the literature concerning multistationarity in this context, mainly from Feliu and Wiuf (J R Soc Interface 10:20130484, 2013), Sadeghimanesh and Feliu (Bull Math Biol 81:2428-2462, 2019), Pérez Millán and Dickenstein (SIAM J Appl Dyn Syst 17(2):1650-1682, 2018), Dickenstein et al. (Bull Math Biol 81:1527-1581, 2019). We establish general results about MESSI systems, which we use to compute the circuits of multistationarity for significant biochemical networks.

A two-strain model of infectious disease spread with asymmetric temporary immunity periods and partial cross-immunity

We introduce a two-strain model with asymmetric temporary immunity periods and partial cross-immunity. We derive explicit conditions for competitive exclusion and coexistence of the strains depending on the strain-specific basic reproduction numbers, temporary immunity periods, and degree of cross-immunity. The results of our bifurcation analysis suggest that, even when two strains share similar basic reproduction numbers and other epidemiological parameters, a disparity in temporary immunity periods and partial or complete cross-immunity can provide a significant competitive advantage. To analyze the dynamics, we introduce a quasi-steady state reduced model which assumes the original strain remains at its endemic steady state. We completely analyze the resulting reduced planar hybrid switching system using linear stability analysis, planar phase-plane analysis, and the Bendixson-Dulac criterion. We validate both the full and reduced models with COVID-19 incidence data, focusing on the Delta (B.1.617.2), Omicron (B.1.1.529), and Kraken (XBB.1.5) variants. These numerical studies suggest that, while early novel strains of COVID-19 had a tendency toward dramatic takeovers and extinction of ancestral strains, more recent strains have the capacity for co-existence.

The unraveling of balanced complexes in metabolic networks

Balanced complexes in biochemical networks are at core of several theoretical and computational approaches that make statements about the properties of the steady states supported by the network. Recent computational approaches have employed balanced complexes to reduce metabolic networks, while ensuring preservation of particular steady-state properties; however, the underlying factors leading to the formation of balanced complexes have not been studied, yet. Here, we present a number of factorizations providing insights in mechanisms that lead to the origins of the corresponding balanced complexes. The proposed factorizations enable us to categorize balanced complexes into four distinct classes, each with specific origins and characteristics. They also provide the means to efficiently determine if a balanced complex in large-scale networks belongs to a particular class from the categorization. The results are obtained under very general conditions and irrespective of the network kinetics, rendering them broadly applicable across variety of network models. Application of the categorization shows that all classes of balanced complexes are present in large-scale metabolic models across all kingdoms of life, therefore paving the way to study their relevance with respect to different properties of steady states supported by these networks.

Multistationarity in Cyclic Sequestration-Transmutation Networks

We consider a natural class of reaction networks which consist of reactions where either two species can inactivate each other (i.e., sequestration), or some species can be transformed into another (i.e., transmutation), in a way that gives rise to a feedback cycle. We completely characterize the capacity of multistationarity of these networks. This is especially interesting because such networks provide simple examples of "atoms of multistationarity", i.e., minimal networks that can give rise to multiple positive steady states.

A structural property for reduction of biochemical networks

Large-scale biochemical models are of increasing sizes due to the consideration of interacting organisms and tissues. Model reduction approaches that preserve the flux phenotypes can simplify the analysis and predictions of steady-state metabolic phenotypes. However, existing approaches either restrict functionality of reduced models or do not lead to significant decreases in the number of modelled metabolites. Here, we introduce an approach for model reduction based on the structural property of balancing of complexes that preserves the steady-state fluxes supported by the network and can be efficiently determined at genome scale. Using two large-scale mass-action kinetic models of Escherichia coli, we show that our approach results in a substantial reduction of 99% of metabolites. Applications to genome-scale metabolic models across kingdoms of life result in up to 55% and 85% reduction in the number of metabolites when arbitrary and mass-action kinetics is assumed, respectively. We also show that predictions of the specific growth rate from the reduced models match those based on the original models. Since steady-state flux phenotypes from the original model are preserved in the reduced, the approach paves the way for analysing other metabolic phenotypes in large-scale biochemical networks.

Stochastic Simulations as a Tool for Assessing Signal Fidelity in Gene Expression in Synthetic Promoter Design

Simple Summary

Synthetic biology is an emerging discipline, offering new perspectives in many industrial fields, from pharma and row-material production to renewable energy. Developing synthetic biology applications is often a lengthy and expensive process with extensive and tedious trial-and-error runs. Computational models can direct the engineering of biological circuits in a computer-aided design setting. By providing a virtual lab environment, in silico models of synthetic circuits can contribute to a quantitative understanding of the underlying molecular pathways before a wet-lab implementation. Here, we illustrate this notion from the point of view of signal fidelity and noise relationship. Noise in gene expression can undermine signal fidelity with implications on the well-functioning of the engineered organisms. For our analysis, we use a specific biological circuit that regulates the gene expression in bacterial inorganic phosphate economy. Applications that use this circuit include those in pollutant detection and wastewater treatment. We provide computational models with different levels of molecular detail as virtual labs. We show that inherent fluctuations in the gene expression machinery can be predicted via stochastic simulations to introduce control in the synthetic promoter design process. Our analysis suggests that noise in the system can be alleviated by strong synthetic promoters with slow unbinding rates. Overall, we provide a recipe for the computer-aided design of synthetic promoter libraries with specific signal to noise characteristics.

Abstract

The design and development of synthetic biology applications in a workflow often involve connecting modular components. Whereas computer-aided design tools are picking up in synthetic biology as in other areas of engineering, the methods for verifying the correct functioning of living technologies are still in their infancy. Especially, fine-tuning for the right promoter strength to match the design specifications is often a lengthy and expensive experimental process. In particular, the relationship between signal fidelity and noise in synthetic promoter design can be a key parameter that can affect the healthy functioning of the engineered organism. To this end, based on our previous work on synthetic promoters for the E. coli PhoBR two-component system, we make a case for using chemical reaction network models for computational verification of various promoter designs before a lab implementation. We provide an analysis of this system with extensive stochastic simulations at a single-cell level to assess the signal fidelity and noise relationship. We then show how quasi-steady-state analysis via ordinary differential equations can be used to navigate between models with different levels of detail. We compare stochastic simulations with our full and reduced models by using various metrics for assessing noise. Our analysis suggests that strong promoters with low unbinding rates can act as control tools for filtering out intrinsic noise in the PhoBR context. Our results confirm that even simpler models can be used to determine promoters with specific signal to noise characteristics.

Dynamics of ERK regulation in the processive limit

We consider a model of extracellular signal-regulated kinase regulation by dual-site phosphorylation and dephosphorylation, which exhibits bistability and oscillations, but loses these properties in the limit in which the mechanisms underlying phosphorylation and dephosphorylation become processive. Our results suggest that anywhere along the way to becoming processive, the model remains bistable and oscillatory. More precisely, in simplified versions of the model, precursors to bistability and oscillations (specifically, multistationarity and Hopf bifurcations, respectively) exist at all "processivity levels". Finally, we investigate whether bistability and oscillations can exist together.

Joining and decomposing reaction networks

In systems and synthetic biology, much research has focused on the behavior and design of single pathways, while, more recently, experimental efforts have focused on how cross-talk (coupling two or more pathways) or inhibiting molecular function (isolating one part of the pathway) affects systems-level behavior. However, the theory for tackling these larger systems in general has lagged behind. Here, we analyze how joining networks (e.g., cross-talk) or decomposing networks (e.g., inhibition or knock-outs) affects three properties that reaction networks may possess-identifiability (recoverability of parameter values from data), steady-state invariants (relationships among species concentrations at steady state, used in model selection), and multistationarity (capacity for multiple steady states, which correspond to multiple cell decisions). Specifically, we prove results that clarify, for a network obtained by joining two smaller networks, how properties of the smaller networks can be inferred from or can imply similar properties of the original network. Our proofs use techniques from computational algebraic geometry, including elimination theory and differential algebra.

Oscillations and bistability in a model of ERK regulation

This work concerns the question of how two important dynamical properties, oscillations and bistability, emerge in an important biological signaling network. Specifically, we consider a model for dual-site phosphorylation and dephosphorylation of extracellular signal-regulated kinase (ERK). We prove that oscillations persist even as the model is greatly simplified (reactions are made irreversible and intermediates are removed). Bistability, however, is much less robust-this property is lost when intermediates are removed or even when all reactions are made irreversible. Moreover, bistability is characterized by the presence of two reversible, catalytic reactions: as other reactions are made irreversible, bistability persists as long as one or both of the specified reactions is preserved. Finally, we investigate the maximum number of steady states, aided by a network's "mixed volume" (a concept from convex geometry). Taken together, our results shed light on the question of how oscillations and bistability emerge from a limiting network of the ERK network-namely, the fully processive dual-site network-which is known to be globally stable and therefore lack both oscillations and bistability. Our proofs are enabled by a Hopf bifurcation criterion due to Yang, analyses of Newton polytopes arising from Hurwitz determinants, and recent characterizations of multistationarity for networks having a steady-state parametrization.

Dynamics of Posttranslational Modification Systems: Recent Progress and Future Directions

Posttranslational modification of proteins is important for signal transduction, and hence significant effort has gone toward understanding how posttranslational modification networks process information. This involves, on the theory side, analyzing the dynamical systems arising from such networks. Which networks are, for instance, bistable? Which networks admit sustained oscillations? Which parameter values enable such behaviors? In this Biophysical Perspective, we highlight recent progress in this area and point out some important future directions. Along the way, we summarize several techniques for analyzing general networks, such as eliminating variables to obtain steady-state parameterizations, and harnessing results on how incorporating intermediates affects dynamics.

Multistationarity and Bistability for Fewnomial Chemical Reaction Networks

Bistability and multistationarity are properties of reaction networks linked to switch-like responses and connected to cell memory and cell decision making. Determining whether and when a network exhibits bistability is a hard and open mathematical problem. One successful strategy consists of analyzing small networks and deducing that some of the properties are preserved upon passage to the full network. Motivated by this, we study chemical reaction networks with few chemical complexes. Under mass action kinetics, the steady states of these networks are described by fewnomial systems, that is polynomial systems having few distinct monomials. Such systems of polynomials are often studied in real algebraic geometry by the use of Gale dual systems. Using this Gale duality, we give precise conditions in terms of the reaction rate constants for the number and stability of the steady states of families of reaction networks with one non-flow reaction.

Regions of multistationarity in cascades of Goldbeter-Koshland loops

We consider cascades of enzymatic Goldbeter-Koshland loops (Goldbeter and Koshland in Proc Natl Acad Sci 78(11):6840-6844, 1981) with any number n of layers, for which there exist two layers involving the same phosphatase. Even if the number of variables and the number of conservation laws grow linearly with n, we find explicit regions in reaction rate constant and total conservation constant space for which the associated mass-action kinetics dynamical system is multistationary. Our computations are based on the theoretical results of our companion paper (Bihan, Dickenstein and Giaroli 2018, preprint: arXiv:1807.05157 ) which are inspired by results in real algebraic geometry by Bihan et al. (SIAM J Appl Algebra Geom, 2018).

Graphical requirements for multistationarity in reaction networks and their verification in BioModels

Thomas' necessary conditions for the existence of multiple steady states in gene networks have been proved by Soulé with high generality for dynamical systems defined by differential equations. When applied to (protein) reaction networks however, those conditions do not provide information since they are trivially satisfied as soon as there is a bimolecular or a reversible reaction. Refined graphical requirements have been proposed to deal with such cases. In this paper, we present for the first time a graph rewriting algorithm for checking the refined conditions given by Soliman, and evaluate its practical performance by applying it systematically to the curated branch of the BioModels repository. This algorithm analyzes all reaction networks (of size up to 430 species) in less than 0.05 second per network, and permits to conclude to the absence of multistationarity in 160 networks over 506. The short computation times obtained in this graphical approach are in sharp contrast to the Jacobian-based symbolic computation approach. We also discuss the case of one extra graphical condition by arc rewiring that allows us to conclude on 20 more networks of this benchmark but with a high computational cost. Finally, we study with some details the case of phosphorylation cycles and MAPK signalling models which show the importance of modelling the intermediate complexations with the enzymes in order to correctly analyze the multistationarity capabilities of such biochemical reaction networks.

Memory functions reveal structural properties of gene regulatory networks

Gene regulatory networks (GRNs) control cellular function and decision making during tissue development and homeostasis. Mathematical tools based on dynamical systems theory are often used to model these networks, but the size and complexity of these models mean that their behaviour is not always intuitive and the underlying mechanisms can be difficult to decipher. For this reason, methods that simplify and aid exploration of complex networks are necessary. To this end we develop a broadly applicable form of the Zwanzig-Mori projection. By first converting a thermodynamic state ensemble model of gene regulation into mass action reactions we derive a general method that produces a set of time evolution equations for a subset of components of a network. The influence of the rest of the network, the bulk, is captured by memory functions that describe how the subnetwork reacts to its own past state via components in the bulk. These memory functions provide probes of near-steady state dynamics, revealing information not easily accessible otherwise. We illustrate the method on a simple cross-repressive transcriptional motif to show that memory functions not only simplify the analysis of the subnetwork but also have a natural interpretation. We then apply the approach to a GRN from the vertebrate neural tube, a well characterised developmental transcriptional network composed of four interacting transcription factors. The memory functions reveal the function of specific links within the neural tube network and identify features of the regulatory structure that specifically increase the robustness of the network to initial conditions. Taken together, the study provides evidence that Zwanzig-Mori projections offer powerful and effective tools for simplifying and exploring the behaviour of GRNs.

Gene regulatory networks are essential for cell fate specification and function. But the recursive links that comprise these networks often make determining their properties and behaviour complicated. Computational models of these networks can also be difficult to decipher. To reduce the complexity of such models we employ a Zwanzig-Mori projection approach. This allows a system of ordinary differential equations, representing a network, to be reduced to an arbitrary subnetwork consisting of part of the initial network, with the rest of the network (bulk) captured by memory functions. These memory functions account for the bulk by describing signals that return to the subnetwork after some time, having passed through the bulk. We show how this approach can be used to simplify analysis and to probe the behaviour of a gene regulatory network. Applying the method to a transcriptional network in the vertebrate neural tube reveals previously unappreciated properties of the network. By taking advantage of the structure of the memory functions we identify interactions within the network that are unnecessary for sustaining correct patterning. Upon further investigation we find that these interactions are important for conferring robustness to variation in initial conditions. Taken together we demonstrate the validity and applicability of the Zwanzig-Mori projection approach to gene regulatory networks.

Core signalling motif displaying multistability through multi-state enzymes

Bistability, and more generally multistability, is a key system dynamics feature enabling decision-making and memory in cells. Deciphering the molecular determinants of multistability is thus crucial for a better understanding of cellular pathways and their (re)engineering in synthetic biology. Here, we show that a key motif found predominantly in eukaryotic signalling systems, namely a futile signalling cycle, can display bistability when featuring a two-state kinase. We provide necessary and sufficient mathematical conditions on the kinetic parameters of this motif that guarantee the existence of multiple steady states. These conditions foster the intuition that bistability arises as a consequence of competition between the two states of the kinase. Extending from this result, we find that increasing the number of kinase states linearly translates into an increase in the number of steady states in the system. These findings reveal, to our knowledge, a new mechanism for the generation of bistability and multistability in cellular signalling systems. Further the futile cycle featuring a two-state kinase is among the smallest bistable signalling motifs. We show that multi-state kinases and the described competition-based motif are part of several natural signalling systems and thereby could enable them to implement complex information processing through multistability. These results indicate that multi-state kinases in signalling systems are readily exploited by natural evolution and could equally be used by synthetic approaches for the generation of multistable information processing systems at the cellular level.

Identifying parameter regions for multistationarity

Mathematical modelling has become an established tool for studying the dynamics of biological systems. Current applications range from building models that reproduce quantitative data to identifying systems with predefined qualitative features, such as switching behaviour, bistability or oscillations. Mathematically, the latter question amounts to identifying parameter values associated with a given qualitative feature. We introduce a procedure to partition the parameter space of a parameterized system of ordinary differential equations into regions for which the system has a unique or multiple equilibria. The procedure is based on the computation of the Brouwer degree, and it creates a multivariate polynomial with parameter depending coefficients. The signs of the coefficients determine parameter regions with and without multistationarity. A particular strength of the procedure is the avoidance of numerical analysis and parameter sampling. The procedure consists of a number of steps. Each of these steps might be addressed algorithmically using various computer programs and available software, or manually. We demonstrate our procedure on several models of gene transcription and cell signalling, and show that in many cases we obtain a complete partitioning of the parameter space with respect to multistationarity.

Intermediates and Generic Convergence to Equilibria

Known graphical conditions for the generic and global convergence to equilibria of the dynamical system arising from a reaction network are shown to be invariant under the so-called successive removal of intermediates, a systematic procedure to simplify the network, making the graphical conditions considerably easier to check.

Chemical Reaction Network Theory elucidates sources of multistability in interferon signaling

Bistability has important implications in signaling pathways, since it indicates a potential cell decision between alternative outcomes. We present two approaches developed in the framework of the Chemical Reaction Network Theory for easy and efficient search of multiple steady state behavior in signaling networks (both with and without mass conservation), and apply them to search for sources of bistability at different levels of the interferon signaling pathway. Different type I interferon subtypes and/or doses are known to elicit differential bioactivities (ranging from antiviral, antiproliferative to immunomodulatory activities). How different signaling outcomes can be generated through the same receptor and activating the same JAK/STAT pathway is still an open question. Here, we detect bistability at the level of early STAT signaling, showing how two different cell outcomes are achieved under or above a threshold in ligand dose or ligand-receptor affinity. This finding could contribute to explain the differential signaling (antiviral vs apoptotic) depending on interferon dose and subtype (α vs β) observed in type I interferons.

Structural simplification of chemical reaction networks in partial steady states

We study the structural simplification of chemical reaction networks with partial steady state semantics assuming that the concentrations of some but not all species are constant. We present a simplification rule that can eliminate intermediate species that are in partial steady state, while preserving the dynamics of all other species. Our simplification rule can be applied to general reaction networks with some but few restrictions on the possible kinetic laws. We can also simplify reaction networks subject to conservation laws. We prove that our simplification rule is correct when applied to a module of a reaction network, as long as the partial steady state is assumed with respect to the complete network. Michaelis-Menten's simplification rule for enzymatic reactions falls out as a special case. We have implemented an algorithm that applies our simplification rules repeatedly and applied it to reaction networks from systems biology.

Intermediates, catalysts, persistence, and boundary steady states

For dynamical systems arising from chemical reaction networks, persistence is the property that each species concentration remains positively bounded away from zero, as long as species concentrations were all positive in the beginning. We describe two graphical procedures for simplifying reaction networks without breaking known necessary or sufficient conditions for persistence, by iteratively removing so-called intermediates and catalysts from the network. The procedures are easy to apply and, in many cases, lead to highly simplified network structures, such as monomolecular networks. For specific classes of reaction networks, we show that these conditions for persistence are equivalent to one another. Furthermore, they can also be characterized by easily checkable strong connectivity properties of a related graph. In particular, this is the case for (conservative) monomolecular networks, as well as cascades of a large class of post-translational modification systems (of which the MAPK cascade and the n-site futile cycle are prominent examples). Since one of the aforementioned sufficient conditions for persistence precludes the existence of boundary steady states, our method also provides a graphical tool to check for that.

Graphical reduction of reaction networks by linear elimination of species

The quasi-steady state approximation and time-scale separation are commonly applied methods to simplify models of biochemical reaction networks based on ordinary differential equations (ODEs). The concentrations of the "fast" species are assumed effectively to be at steady state with respect to the "slow" species. Under this assumption the steady state equations can be used to eliminate the "fast" variables and a new ODE system with only the slow species can be obtained. We interpret a reduced system obtained by time-scale separation as the ODE system arising from a unique reaction network, by identification of a set of reactions and the corresponding rate functions. The procedure is graphically based and can easily be worked out by hand for small networks. For larger networks, we provide a pseudo-algorithm. We study properties of the reduced network, its kinetics and conservation laws, and show that the kinetics of the reduced network fulfil realistic assumptions, provided the original network does. We illustrate our results using biological examples such as substrate mechanisms, post-translational modification systems and networks with intermediates (transient) steps.

Unlimited multistability and Boolean logic in microbial signalling

The ability to map environmental signals onto distinct internal physiological states or programmes is critical for single-celled microbes. A crucial systems dynamics feature underpinning such ability is multistability. While unlimited multistability is known to arise from multi-site phosphorylation seen in the signalling networks of eukaryotic cells, a similarly universal mechanism has not been identified in microbial signalling systems. These systems are generally known as two-component systems comprising histidine kinase (HK) receptors and response regulator proteins engaging in phosphotransfer reactions. We develop a mathematical framework for analysing microbial systems with multi-domain HK receptors known as hybrid and unorthodox HKs. We show that these systems embed a simple core network that exhibits multistability, thereby unveiling a novel biochemical mechanism for multistability. We further prove that sharing of downstream components allows a system with n multi-domain hybrid HKs to attain 3n steady states. We find that such systems, when sensing distinct signals, can readily implement Boolean logic functions on these signals. Using two experimentally studied examples of two-component systems implementing hybrid HKs, we show that bistability and implementation of logic functions are possible under biologically feasible reaction rates. Furthermore, we show that all sequenced microbial genomes contain significant numbers of hybrid and unorthodox HKs, and some genomes have a larger fraction of these proteins compared with regular HKs. Microbial cells are thus theoretically unbounded in mapping distinct environmental signals onto distinct physiological states and perform complex computations on them. These findings facilitate the understanding of natural two-component systems and allow their engineering through synthetic biology.

Morphisms of reaction networks that couple structure to function

BACKGROUND

The mechanisms underlying complex biological systems are routinely represented as networks. Network kinetics is widely studied, and so is the connection between network structure and behavior. However, similarity of mechanism is better revealed by relationships between network structures.

RESULTS

We define morphisms (mappings) between reaction networks that establish structural connections between them. Some morphisms imply kinetic similarity, and yet their properties can be checked statically on the structure of the networks. In particular we can determine statically that a complex network will emulate a simpler network: it will reproduce its kinetics for all corresponding choices of reaction rates and initial conditions. We use this property to relate the kinetics of many common biological networks of different sizes, also relating them to a fundamental population algorithm.

CONCLUSIONS

Structural similarity between reaction networks can be revealed by network morphisms, elucidating mechanistic and functional aspects of complex networks in terms of simpler networks.