Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlögl model revisited

Estimator of entropy production for partially accessible Markov networks based on the observation of blurred transitions

A central task in stochastic thermodynamics is the estimation of entropy production for partially accessible Markov networks. We establish an effective transition-based description for such networks with transitions that are not distinguishable and therefore blurred for an external observer. We demonstrate that, in contrast to a description based on fully resolved transitions, this effective description is typically non-Markovian at any point in time. Starting from an information-theoretic bound, we derive an operationally accessible entropy estimator for this observation scenario. We illustrate the operational relevance and the quality of this entropy estimator with a numerical analysis of various representative examples.

Nonequilibrium fluctuations of chemical reaction networks at criticality: The Schlögl model as paradigmatic case

Chemical reaction networks can undergo nonequilibrium phase transitions upon variation in external control parameters, such as the chemical potential of a species. We investigate the flux in the associated chemostats that is proportional to the entropy production and its critical fluctuations within the Schlögl model. Numerical simulations show that the corresponding diffusion coefficient diverges at the critical point as a function of system size. In the vicinity of the critical point, the diffusion coefficient follows a scaling form. We develop an analytical approach based on the chemical Langevin equation and van Kampen's system size expansion that yields the corresponding exponents in the monostable regime. In the bistable regime, we rely on a two-state approximation in order to analytically describe the critical behavior.

Universal slow dynamics of chemical reaction networks

Understanding the emergent behavior of chemical reaction networks (CRNs) is a fundamental aspect of biology and its origin from inanimate matter. A closed CRN monotonically tends to thermal equilibrium, but when it is opened to external reservoirs, a range of behaviors is possible, including transition to a new equilibrium state, a nonequilibrium state, or indefinite growth. This study shows that slowly driven CRNs are governed by the conserved quantities of the closed system, which are generally far fewer in number than the species. Considering both deterministic and stochastic dynamics, a universal slow-dynamics equation is derived with singular perturbation methods and is shown to be thermodynamically consistent. The slow dynamics is highly robust against microscopic details of the network, which may be unknown in practical situations. In particular, nonequilibrium states of realistic large CRNs can be sought without knowledge of bulk reaction rates. The framework is successfully tested against a suite of networks of increasing complexity and argued to be relevant in the treatment of open CRNs as chemical machines.

Momentum transport of morphological instability in fluid displacement with changes in viscosity

Saffman-Taylor instability exhibits a stepwise unstable morphology from a stable interface to viscous fingering, eventually leading to tip splitting. The nonlinear dynamics of the destabilized interface depends on various flow properties. However, the physicochemical mechanism that determines the transition point of the flow state is unclear. We studied the interfacial instability transition in miscible displacement from a thermodynamic perspective by calculating the momentum transport and entropy production. Using numerical analysis based on Darcy's law coupled with the convection-diffusion equation, the observed flux-dependent flow state transitions were attributed to the selection of the flow state with a higher entropy production.

Macrophage phenotype transitions in a stochastic gene-regulatory network model

Polarization is the process by which a macrophage cell commits to a phenotype based on external signal stimulation. To know how this process is affected by random fluctuations and events within a cell is of utmost importance to better understand the underlying dynamics and predict possible phenotype transitions. For this purpose, we develop a stochastic modeling approach for the macrophage polarization process. We classify phenotype states using the Robust Perron Cluster Analysis and quantify transition pathways and probabilities by applying Transition Path Theory. Depending on the model parameters, we identify four bistable and one tristable phenotype configuration. We find that bistable transitions are fast but their states less robust. In contrast, phenotype transitions in the tristable situation have a comparatively long time duration, which reflects the robustness of the states. The results indicate parallels in the overall transition behavior of macrophage cells with other heterogeneous and plastic cell types, such as cancer cells. Our approach allows for a probabilistic interpretation of macrophage phenotype transitions and biological inference on phenotype robustness. In general, the methodology can easily be adapted to other systems where random state switches are known to occur.

Theoretical investigation of functional responses of bio-molecular assembly networks

Cooperative protein-protein and protein-DNA interactions form programmable complex assemblies, often performing non-linear gene regulatory operations involved in signal transductions and cell fate determination. The apparent structure of those complex assemblies is very similar, but their functional response strongly depends on the topology of the protein-DNA interaction networks. Here, we demonstrate how the coordinated self-assembly creates gene regulatory network motifs that corroborate the existence of a precise functional response at the molecular level using thermodynamic and dynamic analyses. Our theoretical and Monte Carlo simulations show that a complex network of interactions can form a decision-making loop, such as feedback and feed-forward circuits, only by a few molecular mechanisms. We characterize each possible network of interactions by systematic variations of free energy parameters associated with the binding among biomolecules and DNA looping. We also find that the higher-order networks exhibit alternative steady states from the stochastic dynamics of each network. We capture this signature by calculating stochastic potentials and attributing their multi-stability features. We validate our findings against the Gal promoter system in yeast cells. Overall, we show that the network topology is vital in phenotype diversity in regulatory circuits.

Trade-offs between number fluctuations and response in nonequilibrium chemical reaction networks

We study the response of chemical reaction networks driven far from equilibrium to logarithmic perturbations of reaction rates. The response of the mean number of a chemical species is observed to be quantitively limited by number fluctuations and the maximum thermodynamic driving force. We prove these trade-offs for linear chemical reaction networks and a class of nonlinear chemical reaction networks with a single chemical species. Numerical results for several model systems support the conclusion that these trade-offs continue to hold for a broad class of chemical reaction networks, though their precise form appears to sensitively depend on the deficiency of the network.

Foraging Dynamics and Entropy Production in a Simulated Proto-Cell

All organisms depend on a supply of energetic resources to power behavior and the irreversible entropy-producing processes that sustain them. Dissipative structure theory has often been a source of inspiration for better understanding the thermodynamics of biology, yet real organisms are inordinately more complex than most laboratory systems. Here we report on a simulated chemical dissipative structure that operates as a proto cell. The simulated swimmer moves through a 1D environment collecting resources that drive a nonlinear reaction network interior to the swimmer. The model minimally represents properties of a simple organism including rudimentary foraging and chemotaxis and an analog of a metabolism in the nonlinear reaction network. We evaluated how dynamical stability of the foraging dynamics (i.e., swimming and chemotaxis) relates to the rate of entropy production. Results suggested a relationship between dynamical steady states and entropy production that was tuned by the relative coordination of foraging and metabolic processes. Results include evidence in support of and contradicting one formulation of a maximum entropy production principle. We discuss the status of this principle and its relevance to biology.

Operant conditioning of stochastic chemical reaction networks

Adapting one's behavior to environmental conditions and past experience is a key trait of living systems. In the biological world, there is evidence for adaptive behaviors such as learning even in naturally occurring, non-neural, single-celled organisms. In the bioengineered world, advances in synthetic cell engineering and biorobotics have created the possibility of implementing lifelike systems engineered from the bottom up. This will require the development of programmable control circuitry for such biomimetic systems that is capable of realizing such non-trivial and adaptive behavior, including modification of subsequent behavior in response to environmental feedback. To this end, we report the design of novel stochastic chemical reaction networks capable of probabilistic decision-making in response to stimuli. We show that a simple chemical reaction network motif can be tuned to produce arbitrary decision probabilities when choosing between two or more responses to a stimulus signal. We further show that simple feedback mechanisms from the environment can modify these probabilities over time, enabling the system to adapt its behavior dynamically in response to positive or negative reinforcement based on its decisions. This system thus acts as a form of operant conditioning of the chemical circuit, in the sense that feedback provided based on decisions taken by the circuit form the basis of the learning process. Our work thus demonstrates that simple chemical systems can be used to implement lifelike behavior in engineered biomimetic systems.

Efficient Bayesian inference for stochastic agent-based models

The modelling of many real-world problems relies on computationally heavy simulations of randomly interacting individuals or agents. However, the values of the parameters that underlie the interactions between agents are typically poorly known, and hence they need to be inferred from macroscopic observations of the system. Since statistical inference rests on repeated simulations to sample the parameter space, the high computational expense of these simulations can become a stumbling block. In this paper, we compare two ways to mitigate this issue in a Bayesian setting through the use of machine learning methods: One approach is to construct lightweight surrogate models to substitute the simulations used in inference. Alternatively, one might altogether circumvent the need for Bayesian sampling schemes and directly estimate the posterior distribution. We focus on stochastic simulations that track autonomous agents and present two case studies: tumour growths and the spread of infectious diseases. We demonstrate that good accuracy in inference can be achieved with a relatively small number of simulations, making our machine learning approaches orders of magnitude faster than classical simulation-based methods that rely on sampling the parameter space. However, we find that while some methods generally produce more robust results than others, no algorithm offers a one-size-fits-all solution when attempting to infer model parameters from observations. Instead, one must choose the inference technique with the specific real-world application in mind. The stochastic nature of the considered real-world phenomena poses an additional challenge that can become insurmountable for some approaches. Overall, we find machine learning approaches that create direct inference machines to be promising for real-world applications. We present our findings as general guidelines for modelling practitioners.

Network design principle for robust oscillatory behaviors with respect to biological noise

Oscillatory behaviors, which are ubiquitous in transcriptional regulatory networks, are often subject to inevitable biological noise. Thus, a natural question is how transcriptional regulatory networks can robustly achieve accurate oscillation in the presence of biological noise. Here, we search all two- and three-node transcriptional regulatory network topologies for those robustly capable of accurate oscillation against the parameter variability (extrinsic noise) or stochasticity of chemical reactions (intrinsic noise). We find that, no matter what source of the noise is applied, the topologies containing the repressilator with positive autoregulation show higher robustness of accurate oscillation than those containing the activator-inhibitor oscillator, and additional positive autoregulation enhances the robustness against noise. Nevertheless, the attenuation of different sources of noise is governed by distinct mechanisms: the parameter variability is buffered by the long period, while the stochasticity of chemical reactions is filtered by the high amplitude. Furthermore, we analyze the noise of a synthetic human nuclear factor κB (NF-κB) signaling network by varying three different topologies and verify that the addition of a repressilator to the activator-inhibitor oscillator, which leads to the emergence of high-robustness motif—the repressilator with positive autoregulation—improves the oscillation accuracy in comparison to the topology with only an activator-inhibitor oscillator. These design principles may be applicable to other oscillatory circuits.

Current fluctuations in nonequilibrium discontinuous phase transitions

Discontinuous phase transitions out of equilibrium can be characterized by the behavior of macroscopic stochastic currents. But while much is known about the average current, the situation is much less understood for higher statistics. In this paper, we address the consequences of the diverging metastability lifetime-a hallmark of discontinuous transitions-in the fluctuations of arbitrary thermodynamic currents, including the entropy production. In particular, we center our discussion on the conditional statistics, given which phase the system is in. We highlight the interplay between integration window and metastability lifetime, which is not manifested in the average current, but strongly influences the fluctuations. We introduce conditional currents and find, among other predictions, their connection to average and scaled variance through a finite-time version of large deviation theory and a minimal model. Our results are then further verified in two paradigmatic models of discontinuous transitions: Schlögl's model of chemical reactions, and a 12-state Potts model subject to two baths at different temperatures.

Detailed balance, local detailed balance, and global potential for stochastic chemical reaction networks

Detailed balance of a chemical reaction network can be defined in several different ways. Here we investigate the relationship among four types of detailed balance conditions: deterministic, stochastic, local, and zero-order local detailed balance. We show that the four types of detailed balance are equivalent when different reactions lead to different species changes and are not equivalent when some different reactions lead to the same species change. Under the condition of local detailed balance, we further show that the system has a global potential defined over the whole space, which plays a central role in the large deviation theory and the Freidlin–Wentzell-type metastability theory of chemical reaction networks. Finally, we provide a new sufficient condition for stochastic detailed balance, which is applied to construct a class of high-dimensional chemical reaction networks that both satisfies stochastic detailed balance and displays multistability.

Noise distorts the epigenetic landscape and shapes cell-fate decisions

The Waddington epigenetic landscape has become an iconic representation of the cellular differentiation process. Recent single-cell transcriptomic data provide new opportunities for quantifying this originally conceptual tool, offering insight into the gene regulatory networks underlying cellular development. While many methods for constructing the landscape have been proposed, by far the most commonly employed approach is based on computing the landscape as the negative logarithm of the steady-state probability distribution. Here, we use simple models to highlight the complexities and limitations that arise when reconstructing the potential landscape in the presence of stochastic fluctuations. We consider how the landscape changes in accordance with different stochastic systems and show that it is the subtle interplay between the deterministic and stochastic components of the system that ultimately shapes the landscape. We further discuss how the presence of noise has important implications for the identifiability of the regulatory dynamics from experimental data. A record of this paper's transparent peer review process is included in the supplemental information.

Exact Probability Landscapes of Stochastic Phenotype Switching in Feed-Forward Loops: Phase Diagrams of Multimodality

Feed-forward loops (FFLs) are among the most ubiquitously found motifs of reaction networks in nature. However, little is known about their stochastic behavior and the variety of network phenotypes they can exhibit. In this study, we provide full characterizations of the properties of stochastic multimodality of FFLs, and how switching between different network phenotypes are controlled. We have computed the exact steady-state probability landscapes of all eight types of coherent and incoherent FFLs using the finite-butter Accurate Chemical Master Equation (ACME) algorithm, and quantified the exact topological features of their high-dimensional probability landscapes using persistent homology. Through analysis of the degree of multimodality for each of a set of 10,812 probability landscapes, where each landscape resides over 10⁵–10⁶ microstates, we have constructed comprehensive phase diagrams of all relevant behavior of FFL multimodality over broad ranges of input and regulation intensities, as well as different regimes of promoter binding dynamics. In addition, we have quantified the topological sensitivity of the multimodality of the landscapes to regulation intensities. Our results show that with slow binding and unbinding dynamics of transcription factor to promoter, FFLs exhibit strong stochastic behavior that is very different from what would be inferred from deterministic models. In addition, input intensity play major roles in the phenotypes of FFLs: At weak input intensity, FFL exhibit monomodality, but strong input intensity may result in up to 6 stable phenotypes. Furthermore, we found that gene duplication can enlarge stable regions of specific multimodalities and enrich the phenotypic diversity of FFL networks, providing means for cells toward better adaptation to changing environment. Our results are directly applicable to analysis of behavior of FFLs in biological processes such as stem cell differentiation and for design of synthetic networks when certain phenotypic behavior is desired.

Improved estimations of stochastic chemical kinetics by finite-state expansion

Stochastic reaction networks are a fundamental model to describe interactions between species where random fluctuations are relevant. The master equation provides the evolution of the probability distribution across the discrete state space consisting of vectors of population counts for each species. However, since its exact solution is often elusive, several analytical approximations have been proposed. The deterministic rate equation (DRE) gives a macroscopic approximation as a compact system of differential equations that estimate the average populations for each species, but it may be inaccurate in the case of nonlinear interaction dynamics. Here we propose finite-state expansion (FSE), an analytical method mediating between the microscopic and the macroscopic interpretations of a stochastic reaction network by coupling the master equation dynamics of a chosen subset of the discrete state space with the mean population dynamics of the DRE. An algorithm translates a network into an expanded one where each discrete state is represented as a further distinct species. This translation exactly preserves the stochastic dynamics, but the DRE of the expanded network can be interpreted as a correction to the original one. The effectiveness of FSE is demonstrated in models that challenge state-of-the-art techniques due to intrinsic noise, multi-scale populations and multi-stability.

Direct statistical inference for finite Markov jump processes via the matrix exponential

Given noisy, partial observations of a time-homogeneous, finite-statespace Markov chain, conceptually simple, direct statistical inference is available, in theory, via its rate matrix, or infinitesimal generator, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathsf {Q}}$$\end{document}Q, since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\exp ({\mathsf {Q}}t)$$\end{document}exp(Qt) is the transition matrix over time t. However, perhaps because of inadequate tools for matrix exponentiation in programming languages commonly used amongst statisticians or a belief that the necessary calculations are prohibitively expensive, statistical inference for continuous-time Markov chains with a large but finite state space is typically conducted via particle MCMC or other relatively complex inference schemes. When, as in many applications \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathsf {Q}}$$\end{document}Q arises from a reaction network, it is usually sparse. We describe variations on known algorithms which allow fast, robust and accurate evaluation of the product of a non-negative vector with the exponential of a large, sparse rate matrix. Our implementation uses relatively recently developed, efficient, linear algebra tools that take advantage of such sparsity. We demonstrate the straightforward statistical application of the key algorithm on a model for the mixing of two alleles in a population and on the Susceptible-Infectious-Removed epidemic model.

Entropic analysis of bistability and the general evolution criterion

An open flow bistable Schlögl model in a well-stirred isothermal reaction volume V showing the balance of the entropy production and the exchange entropy on the non-equilibrium stationary states.

Fluctuating-rate model with multiple gene states

Multiple phenotypic states of single cells often co-exist in the presence of positive feedbacks. Stochastic gene-state switchings and low copy numbers of proteins in single cells cause considerable fluctuations. The chemical master equation (CME) is a powerful tool that describes the dynamics of single cells, but it may be overly complicated. Among many simplified models, a fluctuating-rate (FR) model has been proposed recently to approximate the full CME model in the realistic intermediate region of gene-state switchings. However, only the scenario with two gene states has been carefully analysed. In this paper, we generalise the FR model to the case with multiple gene states, in which the mathematical derivation becomes more complicated. The leading order of fluctuations around each phenotypic state, as well as the transition rates between phenotypic states, in the intermediate gene-state switching region is characterized by the rate function of the stationary distribution of the FR model in the Freidlin-Wentzell-type large deviation principle (LDP). Under certain reasonable assumptions, we show that the derivative of the rate function is equal to the unique nontrivial solution of a dominant generalised eigenvalue problem, leading to a new numerical algorithm for obtaining the LDP rate function directly. Furthermore, we prove the Lyapunov property of the rate function for the corresponding deterministic mean-field dynamics. Finally, through a tristable example, we show that the local fluctuations (the asymptotic variance of the stationary distribution at each phenotypic state) in the intermediate and rapid regions of gene-state switchings are different. Finally, a tri-stable example is constructed to illustrate the validity of our theory.

Optimizing Our Patients' Entropy Production as Therapy? Hypotheses Originating from the Physics of Physiology

Understanding how nature drives entropy production offers novel insights regarding patient care. Whilst energy is always preserved and energy gradients irreversibly dissipate (thus producing entropy), increasing evidence suggests that they do so in the most optimal means possible. For living complex non-equilibrium systems to create a healthy internal emergent order, they must continuously produce entropy over time. The Maximum Entropy Production Principle (MEPP) highlights nature's drive for non-equilibrium systems to augment their entropy production if possible. This physical drive is hypothesized to be responsible for the spontaneous formation of fractal structures in space (e.g., multi-scale self-similar tree-like vascular structures that optimize delivery to and clearance from an organ system) and time (e.g., complex heart and respiratory rate variability); both are ubiquitous and essential for physiology and health. Second, human entropy production, measured by heat production divided by temperature, is hypothesized to relate to both metabolism and consciousness, dissipating oxidative energy gradients and reducing information into meaning and memory, respectively. Third, both MEPP and natural selection are hypothesized to drive enhanced functioning and adaptability, selecting states with robust basilar entropy production, as well as the capacity to enhance entropy production in response to exercise, heat stress, and illness. Finally, a targeted focus on optimizing our patients' entropy production has the potential to improve health and clinical outcomes. With the implications of developing a novel understanding of health, illness, and treatment strategies, further exploration of this uncharted ground will offer value.

Exponential volume dependence of entropy-current fluctuations at first-order phase transitions in chemical reaction networks

In chemical reaction networks, bistability can only occur far from equilibrium. It is associated with a first-order phase transition where the control parameter is the thermodynamic force. At the bistable point, the entropy production is known to be discontinuous with respect to the thermodynamic force. We show that the fluctuations of the entropy production have an exponential volume-dependence when the system is bistable. At the phase transition, the exponential prefactor is the height of the effective potential barrier between the two fixed-points. Our results obtained for Schlögl's model can be extended to any chemical network.

Cell-to-Cell Information at a Feedback-Induced Bifurcation Point

A ubiquitous way that cells share information is by exchanging molecules. Yet, the fundamental ways that this information exchange is influenced by intracellular dynamics remain unclear. Here we use information theory to investigate a simple model of two interacting cells with internal feedback. We show that cell-to-cell molecule exchange induces a collective two-cell critical point and that the mutual information between the cells peaks at this critical point. Information can remain large far from the critical point on a manifold of cellular states but scales logarithmically with the correlation time of the system, resulting in an information-correlation time trade-off. This trade-off is strictly imposed, suggesting the correlation time as a proxy for the mutual information.

A practical guide to pseudo-marginal methods for computational inference in systems biology

For many stochastic models of interest in systems biology, such as those describing biochemical reaction networks, exact quantification of parameter uncertainty through statistical inference is intractable. Likelihood-free computational inference techniques enable parameter inference when the likelihood function for the model is intractable but the generation of many sample paths is feasible through stochastic simulation of the forward problem. The most common likelihood-free method in systems biology is approximate Bayesian computation that accepts parameters that result in low discrepancy between stochastic simulations and measured data. However, it can be difficult to assess how the accuracy of the resulting inferences are affected by the choice of acceptance threshold and discrepancy function. The pseudo-marginal approach is an alternative likelihood-free inference method that utilises a Monte Carlo estimate of the likelihood function. This approach has several advantages, particularly in the context of noisy, partially observed, time-course data typical in biochemical reaction network studies. Specifically, the pseudo-marginal approach facilitates exact inference and uncertainty quantification, and may be efficiently combined with particle filters for low variance, high-accuracy likelihood estimation. In this review, we provide a practical introduction to the pseudo-marginal approach using inference for biochemical reaction networks as a series of case studies. Implementations of key algorithms and examples are provided using the Julia programming language; a high performance, open source programming language for scientific computing (https://github.com/davidwarne/Warne2019_GuideToPseudoMarginal).

Optimal control of protein copy number

Cell-cell communication is often achieved by secreted signaling molecules that bind membrane-bound receptors. A common class of such receptors are G-protein coupled receptors, where extracellular binding induces changes on the membrane affinity near the receptor for certain cytosolic proteins, effectively altering their chemical potential. We analyze the minimum-dissipation schedules for dynamically changing chemical potential to induce steady-state changes in protein copy-number distributions, and illustrate with analytic solutions for linear chemical reaction networks. Protocols that change chemical potential on biologically relevant timescales are experimentally accessible using optogenetic manipulations, and our framework provides nontrivial predictions about functional dynamical cell-cell interactions.

Thermodynamics of switching in multistable non-equilibrium systems

Multistable non-equilibrium systems are abundant outcomes of nonlinear dynamics with feedback, but still relatively little is known about what determines the stability of the steady states and their switching rates in terms of entropy and entropy production. Here, we will link fluctuation theorems for the entropy production along trajectories with the action obtainable from the Freidlin-Wentzell theorem to elucidate the thermodynamics of switching between states in the large volume limit of multistable systems. We find that the entropy production at steady state plays no role, but the entropy production during switching is key. Steady-state entropy and diffusive noise strength can be neglected in this limit. The relevance to biological, ecological, and climate models is apparent.

Finite state projection for approximating the stationary solution to the chemical master equation using reaction rate equations

When modeling a physical system using a Markov chain, it is often instructive to compute its probability distribution at statistical equilibrium, thereby gaining insight into the stationary, or long-term, behavior of the system. Computing such a distribution directly is problematic when the state space of the system is large. Here, we look at the case of a chemical reaction system that models the dynamics of cellular processes, where it has become popular to constrain the computational burden by using a finite state projection, which aims only to capture the most likely states of the system, rather than every possible state. We propose an efficient method to further narrow these states to those that remain highly probable in the long run, after the transient behavior of the system has dissipated. Our strategy is to quickly estimate the local maxima of the stationary distribution using the reaction rate formulation, which is of considerably smaller size than the full-blown chemical master equation, and from there develop adaptive schemes to profile the distribution around the maxima. We include numerical tests that show the efficiency of our approach.

Interlinked GTPase cascades provide a motif for both robust switches and oscillators

GTPases regulate a wide range of cellular processes, such as intracellular vesicular transport, signal transduction and protein translation. These hydrolase enzymes operate as biochemical switches by toggling between an active guanosine triphosphate (GTP)-bound state and an inactive guanosine diphosphate (GDP)-bound state. We compare two network motifs, a single-species switch and an interlinked cascade that consists of two species coupled through positive and negative feedback loops. We find that interlinked cascades are closer to the ideal all-or-none switch and are more robust against fluctuating signals. While the single-species switch can only achieve bistability, interlinked cascades can be converted into oscillators by tuning the cofactor concentrations, which catalyse the activity of the cascade. These regimes can only be achieved with sufficient chemical driving provided by GTP hydrolysis. In this study, we present a thermodynamically consistent model that can achieve bistability and oscillations with the same feedback motif.

Critical slowing down in biochemical networks with feedback

Near a bifurcation point, the response time of a system is expected to diverge due to the phenomenon of critical slowing down. We investigate critical slowing down in well-mixed stochastic models of biochemical feedback by exploiting a mapping to the mean-field Ising universality class. We analyze the responses to a sudden quench and to continuous driving in the model parameters. In the latter case, we demonstrate that our class of models exhibits the Kibble-Zurek collapse, which predicts the scaling of hysteresis in cellular responses to gradual perturbations. We discuss the implications of our results in terms of the tradeoff between a precise and a fast response. Finally, we use our mapping to quantify critical slowing down in T cells, where the addition of a drug is equivalent to a sudden quench in parameter space.

Thermodynamic analysis of thermal convection based on entropy production

Flow patterns have a tendency to break the symmetry of an initial state of a system and form another spatiotemporal pattern when the system is driven far from equilibrium by temperature difference. For an annular channel, the axially symmetric flow becomes unstable beyond a given temperature difference threshold imposed in the system, leading to rotational oscillating waves. Many researchers have investigated this transition via linear stability analysis using the fundamental conservation equations and the generic model amplitude equation, i.e., the complex Ginzburg-Landau equation. Here, we present a quantitative study conducted of the thermal convection transition using thermodynamic analysis based on the maximum entropy production principle. Our analysis results reveal that the fluid system under nonequilibrium maximizes the entropy production induced by the thermodynamic flux in a direction perpendicular to the temperature difference. Further, we show that the thermodynamic flux as well as the entropy production can uniquely specify the thermodynamic states of the entire fluid system and propose an entropy production selection rule that can be used to specify the thermodynamic state of a nonequilibrium system.

Universality of biochemical feedback and its application to immune cells

We map a class of well-mixed stochastic models of biochemical feedback in steady state to the mean-field Ising model near the critical point. The mapping provides an effective temperature, magnetic field, order parameter, and heat capacity that can be extracted from biological data without fitting or knowledge of the underlying molecular details. We demonstrate this procedure on fluorescence data from mouse T cells, which reveals distinctions between how the cells respond to different drugs. We also show that the heat capacity allows inference of the absolute molecule number from fluorescence intensity. We explain this result in terms of the underlying fluctuations, and we demonstrate the generality of our work.

A comparison of deterministic and stochastic approaches for sensitivity analysis in computational systems biology

With the recent rising application of mathematical models in the field of computational systems biology, the interest in sensitivity analysis methods had increased. The stochastic approach, based on chemical master equations, and the deterministic approach, based on ordinary differential equations (ODEs), are the two main approaches for analyzing mathematical models of biochemical systems. In this work, the performance of these approaches to compute sensitivity coefficients is explored in situations where stochastic and deterministic simulation can potentially provide different results (systems with unstable steady states, oscillators with population extinction and bistable systems). We consider two methods in the deterministic approach, namely the direct differential method and the finite difference method, and five methods in the stochastic approach, namely the Girsanov transformation, the independent random number method, the common random number method, the coupled finite difference method and the rejection-based finite difference method. The reviewed methods are compared in terms of sensitivity values and computational time to identify differences in outcome that can highlight conditions in which one approach performs better than the other.

Discrete flux and velocity fields of probability and their global maps in reaction systems

Stochasticity plays important roles in reaction systems. Vector fields of probability flux and velocity characterize time-varying and steady-state properties of these systems, including high probability paths, barriers, checkpoints among different stable regions, as well as mechanisms of dynamic switching among them. However, conventional fluxes on continuous space are ill-defined and are problematic when at the boundaries of the state space or when copy numbers are small. By re-defining the derivative and divergence operators based on the discrete nature of reactions, we introduce new formulations of discrete fluxes. Our flux model fully accounts for the discreetness of both the state space and the jump processes of reactions. The reactional discrete flux satisfies the continuity equation and describes the behavior of the system evolving along directions of reactions. The species discrete flux directly describes the dynamic behavior in the state space of the reactants such as the transfer of probability mass. With the relationship between these two fluxes specified, we show how to construct time-evolving and steady-state global flow-maps of probability flux and velocity in the directions of every species at every microstate and how they are related to the outflow and inflow of probability fluxes when tracing out reaction trajectories. We also describe how to impose proper conditions enabling exact quantification of flux and velocity in the boundary regions, without the difficulty of enforcing artificial reflecting conditions. We illustrate the computation of probability flux and velocity using three model systems, namely, the birth-death process, the bistable Schlögl model, and the oscillating Schnakenberg model.

Non-Equilibrium Thermodynamics and Stochastic Dynamics of a Bistable Catalytic Surface Reaction

Catalytic surface reaction networks exhibit nonlinear dissipative phenomena, such as bistability. Macroscopic rate law descriptions predict that the reaction system resides on one of the two steady-state branches of the bistable region for an indefinite period of time. However, the smaller the catalytic surface, the greater the influence of coverage fluctuations, given that their amplitude normally scales as the square root of the system size. Thus, one can observe fluctuation-induced transitions between the steady-states. In this work, a model for the bistable catalytic CO oxidation on small surfaces is studied. After a brief introduction of the average stochastic modelling framework and its corresponding deterministic limit, we discuss the non-equilibrium conditions necessary for bistability. The entropy production rate, an important thermodynamic quantity measuring dissipation in a system, is compared across the two approaches. We conclude that, in our catalytic model, the most favorable non-equilibrium steady state is not necessary the state with the maximum or minimum entropy production rate.

On Differences between Deterministic and Stochastic Models of Chemical Reactions: Schlögl Solved with ZI-Closure

Deterministic and stochastic models of chemical reaction kinetics can give starkly different results when the deterministic model exhibits more than one stable solution. For example, in the stochastic Schlögl model, the bimodal stationary probability distribution collapses to a unimodal distribution when the system size increases, even for kinetic constant values that result in two distinct stable solutions in the deterministic Schlögl model. Using zero-information (ZI) closure scheme, an algorithm for solving chemical master equations, we compute stationary probability distributions for varying system sizes of the Schlögl model. With ZI-closure, system sizes can be studied that have been previously unattainable by stochastic simulation algorithms. We observe and quantify paradoxical discrepancies between stochastic and deterministic models and explain this behavior by postulating that the entropy of non-equilibrium steady states (NESS) is maximum.

A landscape view on the interplay between EMT and cancer metastasis

The epithelial-mesenchymal transition (EMT) is a basic developmental process that converts epithelial cells to mesenchymal cells. Although EMT might promote cancer metastasis, the molecular mechanisms for it remain to be fully clarified. To address this issue, we constructed an EMT-metastasis gene regulatory network model and quantified the potential landscape of cancer metastasis-promoting system computationally. We identified four steady-state attractors on the landscape, which separately characterize anti-metastatic (A), metastatic (M), and two other intermediate (I1 and I2) cell states. The tetrastable landscape and the existence of intermediate states are consistent with recent single-cell measurements. We identified one of the two intermediate states I1 as the EMT state. From a MAP approach, we found that for metastatic progression cells need to first undergo EMT (enter the I1 state), and then become metastatic (switch from the I1 state to the M state). Specifically, for metastatic progression, EMT genes (such as ZEB) should be activated before metastasis genes (such as BACH1). This suggests that temporal order is important for the activation of cellular programs in biological systems, and provides a possible mechanism of EMT-promoting cancer metastasis. To identify possible therapeutic targets from this landscape view, we performed sensitivity analysis for individual molecular factors, and identified optimal interventions for landscape control. We found that minimizing transition actions more effectively identifies optimal combinations of targets that induce transitions between attractors than single-factor sensitivity analysis. Overall, the landscape view not only suggests that intermediate states increase plasticity during cell fate decisions, providing a possible source for tumor heterogeneity that is critically important in metastatic progress, but also provides a way to identify therapeutic targets for preventing cancer progression.

Phase transition in thermodynamically consistent biochemical oscillators

Biochemical oscillations are ubiquitous in living organisms. In an autonomous system, not influenced by an external signal, they can only occur out of equilibrium. We show that they emerge through a generic nonequilibrium phase transition, with a characteristic qualitative behavior at criticality. The control parameter is the thermodynamic force which must be above a certain threshold for the onset of biochemical oscillations. This critical behavior is characterized by the thermodynamic flux associated with the thermodynamic force, its diffusion coefficient, and the stationary distribution of the oscillating chemical species. We discuss metrics for the precision of biochemical oscillations by comparing two observables, the Fano factor associated with the thermodynamic flux and the number of coherent oscillations. Since the Fano factor can be small even when there are no biochemical oscillations, we argue that the number of coherent oscillations is more appropriate to quantify the precision of biochemical oscillations. Our results are obtained with three thermodynamically consistent versions of known models: the Brusselator, the activator-inhibitor model, and a model for KaiC oscillations.

A mechanism for bistability in glycosylation

Glycosyltransferases are a class of enzymes that catalyse the posttranslational modification of proteins to produce a large number of glycoconjugate acceptors from a limited number of nucleotide-sugar donors. The products of one glycosyltransferase can be the substrates of several other enzymes, causing a combinatorial explosion in the number of possible glycan products. The kinetic behaviour of systems where multiple acceptor substrates compete for a single enzyme is presented, and the case in which high concentrations of an acceptor substrate are inhibitory as a result of abortive complex formation, is shown to result in non-Michaelian kinetics that can lead to bistability in an open system. A kinetic mechanism is proposed that is consistent with the available experimental evidence and provides a possible explanation for conflicting observations on the β-1,4-galactosyltransferases. Abrupt switching between steady states in networks of glycosyltransferase-catalysed reactions may account for the observed changes in glycosyl-epitopes in cancer cells.

Stochastic thermodynamics and entropy production of chemical reaction systems

We investigate the nonequilibrium stationary states of systems consisting of chemical reactions among molecules of several chemical species. To this end, we introduce and develop a stochastic formulation of nonequilibrium thermodynamics of chemical reaction systems based on a master equation defined on the space of microscopic chemical states and on appropriate definitions of entropy and entropy production. The system is in contact with a heat reservoir and is placed out of equilibrium by the contact with particle reservoirs. In our approach, the fluxes of various types, such as the heat and particle fluxes, play a fundamental role in characterizing the nonequilibrium chemical state. We show that the rate of entropy production in the stationary nonequilibrium state is a bilinear form in the affinities and the fluxes of reaction, which are expressed in terms of rate constants and transition rates, respectively. We also show how the description in terms of microscopic states can be reduced to a description in terms of the numbers of particles of each species, from which follows the chemical master equation. As an example, we calculate the rate of entropy production of the first and second Schlögl reaction models.

Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach

We develop numerical methods for stochastic reaction-diffusion systems based on approaches used for fluctuatinghydrodynamics (FHD). For hydrodynamicsystems, the FHD formulation is formally described by stochastic partial differential equations (SPDEs). In the reaction-diffusion systems we consider, our model becomes similar to the reaction-diffusion master equation (RDME) description when our SPDEs are spatially discretized and reactions are modeled as a source term having Poissonfluctuations. However, unlike the RDME, which becomes prohibitively expensive for an increasing number of molecules, our FHD-based description naturally extends from the regime where fluctuations are strong, i.e., each mesoscopic cell has few (reactive) molecules, to regimes with moderate or weak fluctuations, and ultimately to the deterministic limit. By treating diffusion implicitly, we avoid the severe restriction on time step size that limits all methods based on explicit treatments of diffusion and construct numerical methods that are more efficient than RDME methods, without compromising accuracy. Guided by an analysis of the accuracy of the distribution of steady-state fluctuations for the linearized reaction-diffusion model, we construct several two-stage (predictor-corrector) schemes, where diffusion is treated using a stochastic Crank-Nicolson method, and reactions are handled by the stochastic simulation algorithm of Gillespie or a weakly second-order tau leaping method. We find that an implicit midpoint tau leaping scheme attains second-order weak accuracy in the linearized setting and gives an accurate and stable structure factor for a time step size of an order of magnitude larger than the hopping time scale of diffusing molecules. We study the numerical accuracy of our methods for the Schlögl reaction-diffusion model both in and out of thermodynamic equilibrium. We demonstrate and quantify the importance of thermodynamicfluctuations to the formation of a two-dimensional Turing-like pattern and examine the effect of fluctuations on three-dimensional chemical front propagation. By comparing stochastic simulations to deterministic reaction-diffusion simulations, we show that fluctuations accelerate pattern formation in spatially homogeneous systems and lead to a qualitatively different disordered pattern behind a traveling wave.

Processes on the emergent landscapes of biochemical reaction networks and heterogeneous cell population dynamics: differentiation in living matters

The notion of an attractor has been widely employed in thinking about the nonlinear dynamics of organisms and biological phenomena as systems and as processes. The notion of a landscape with valleys and mountains encoding multiple attractors, however, has a rigorous foundation only for closed, thermodynamically non-driven, chemical systems, such as a protein. Recent advances in the theory of nonlinear stochastic dynamical systems and its applications to mesoscopic reaction networks, one reaction at a time, have provided a new basis for a landscape of open, driven biochemical reaction systems under sustained chemostat. The theory is equally applicable not only to intracellular dynamics of biochemical regulatory networks within an individual cell but also to tissue dynamics of heterogeneous interacting cell populations. The landscape for an individual cell, applicable to a population of isogenic non-interacting cells under the same environmental conditions, is defined on the counting space of intracellular chemical compositions x = (x₁,x₂, … ,x_N ) in a cell, where x_ℓ is the concentration of the ℓth biochemical species. Equivalently, for heterogeneous cell population dynamics x_ℓ is the number density of cells of the ℓth cell type. One of the insights derived from the landscape perspective is that the life history of an individual organism, which occurs on the hillsides of a landscape, is nearly deterministic and 'programmed', while population-wise an asynchronous non-equilibrium steady state resides mostly in the lowlands of the landscape. We argue that a dynamic 'blue-sky' bifurcation, as a representation of Waddington's landscape, is a more robust mechanism for a cell fate decision and subsequent differentiation than the widely pictured pitch-fork bifurcation. We revisit, in terms of the chemostatic driving forces upon active, living matter, the notions of near-equilibrium thermodynamic branches versus far-from-equilibrium states. The emergent landscape perspective permits a quantitative discussion of a wide range of biological phenomena as nonlinear, stochastic dynamics.

Entrance and escape dynamics for the typical set

According to the asymptotic equipartition property, sufficiently long sequences of random variables converge to a set that is typical. While the size and probability of this set are central to information theory and statistical mechanics, they can often only be estimated accurately in the asymptotic limit due to the exponential growth in possible sequences. Here we derive a time-inhomogeneous dynamics that constructs the properties of the typical set for all finite length sequences of independent and identically distributed random variables. These dynamics link the finite properties of the typical set to asymptotic results and allow the typical set to be applied to small and transient systems. The main result is a geometric mapping-the triangle map-relating sequences of random variables of length n to those of length n+1. We show that the number of points in this map needed to quantify the properties of the typical set grows linearly with sequence length, despite the exponential growth in the number of typical sequences. We illustrate the framework for the Bernoulli process and the Schlögl model for autocatalytic chemical reactions and demonstrate both the convergence to asymptotic limits and the ability to reproduce exact calculations.

Parallel replica dynamics method for bistable stochastic reaction networks: Simulation and sensitivity analysis

Stochastic reaction networks that exhibit bistable behavior are common in systems biology, materials science, and catalysis. Sampling of stationary distributions is crucial for understanding and characterizing the long-time dynamics of bistable stochastic dynamical systems. However, simulations are often hindered by the insufficient sampling of rare transitions between the two metastable regions. In this paper, we apply the parallel replica method for a continuous time Markov chain in order to improve sampling of the stationary distribution in bistable stochastic reaction networks. The proposed method uses parallel computing to accelerate the sampling of rare transitions. Furthermore, it can be combined with the path-space information bounds for parametric sensitivity analysis. With the proposed methodology, we study three bistable biological networks: the Schlögl model, the genetic switch network, and the enzymatic futile cycle network. We demonstrate the algorithmic speedup achieved in these numerical benchmarks. More significant acceleration is expected when multi-core or graphics processing unit computer architectures and programming tools such as CUDA are employed.

Entropy production selects nonequilibrium states in multistable systems

Far-from-equilibrium thermodynamics underpins the emergence of life, but how has been a long-outstanding puzzle. Best candidate theories based on the maximum entropy production principle could not be unequivocally proven, in part due to complicated physics, unintuitive stochastic thermodynamics, and the existence of alternative theories such as the minimum entropy production principle. Here, we use a simple, analytically solvable, one-dimensional bistable chemical system to demonstrate the validity of the maximum entropy production principle. To generalize to multistable stochastic system, we use the stochastic least-action principle to derive the entropy production and its role in the stability of nonequilibrium steady states. This shows that in a multistable system, all else being equal, the steady state with the highest entropy production is favored, with a number of implications for the evolution of biological, physical, and geological systems.

Energy conservation and maximal entropy production in enzyme reactions

A procedure for maximization of the density of entropy production in a single stationary two-step enzyme reaction is developed. Under the constraints of mass conservation, fixed equilibrium constant of a reaction and fixed products of forward and backward enzyme rate constants the existence of maximum in the density of entropy production is demonstrated. In the state with maximal density of entropy production the optimal enzyme rate constants, the stationary concentrations of the substrate and the product, the stationary product yield as well as the stationary reaction flux are calculated. The test, whether these calculated values of the reaction parameters are consistent with their corresponding measured values, is performed for the enzyme Glucose Isomerase. It is found that calculated and measured rate constants agree within an order of magnitude, whereas the calculated reaction flux and the product yield differ from their corresponding measured values for less than 20 % and 5 %, respectively. This indicates that the enzyme Glucose Isomerase, considered in a non-equilibrium stationary state, as found in experiments using the continuous stirred tank reactors, possibly operates close to the state with the maximum in the density of entropy production.

A minimal model of burst-noise induced bistability

We investigate the influence of intrinsic noise on stable states of a one-dimensional dynamical system that shows in its deterministic version a saddle-node bifurcation between monostable and bistable behaviour. The system is a modified version of the Schlögl model, which is a chemical reaction system with only one type of molecule. The strength of the intrinsic noise is varied without changing the deterministic description by introducing bursts in the autocatalytic production step. We study the transitions between monostable and bistable behavior in this system by evaluating the number of maxima of the stationary probability distribution. We find that changing the size of bursts can destroy and even induce saddle-node bifurcations. This means that a bursty production of molecules can qualitatively change the dynamics of a chemical reaction system even when the deterministic description remains unchanged.