Simulation of biochemical reactions with time-dependent rates by the rejection-based algorithm

Investigating seasonal disease emergence and extinction in stochastic epidemic models

Seasonal disease outbreaks are common in many infectious diseases such as seasonal influenza, Zika, dengue fever, Lyme disease, malaria, and cholera. Seasonal outbreaks are often due to weather patterns affecting pathogens or disease-carrying vectors or by social behavior. We investigate disease emergence and extinction in seasonal stochastic epidemic models. Specifically, we study disease emergence through seasonally varying parameters for transmission, recovery, and vector births and deaths in time-nonhomogeneous Markov chains for SIR, SEIR, and vector-host systems. A branching process approximation of the Markov chain is used to estimate the seasonal probabilities of disease extinction. Several disease outcome measures are used to compare the dynamics in seasonal and constant environments. Numerical investigations illustrate and confirm previous results derived from stochastic epidemic models. Seasonal environments often result in lower probabilities of disease emergence and smaller values of the basic reproduction number than in constant environments, and the time of peak emergence generally precedes the peak time of the seasonal driver. We identify some new results when both transmission and recovery vary seasonally. If the relative amplitude of the recovery exceeds that of transmission or if the periodicity is not synchronized in time, lower average probabilities of disease emergence occur in a constant environment than in a seasonal environment. We also investigate the timing of vector control. This investigation provides new methods and outcome measures to study seasonal infectious disease dynamics and offers new insights into the timing of prevention and control.

Quantifying infectious disease epidemic risks: A practical approach for seasonal pathogens

For many infectious diseases, the risk of outbreaks varies seasonally. If a pathogen is usually absent from a host population, a key public health policy question is whether the pathogen's arrival will initiate local transmission, which depends on the season in which arrival occurs. This question can be addressed by estimating the "probability of a major outbreak" (the probability that introduced cases will initiate sustained local transmission). A standard approach for inferring this probability exists for seasonal pathogens (involving calculating the Case Epidemic Risk; CER) based on the mathematical theory of branching processes. Under that theory, the probability of pathogen extinction is estimated, neglecting depletion of susceptible individuals. The CER is then one minus the extinction probability. However, as we show, if transmission cannot occur for long periods of the year (e.g., over winter or over summer), the pathogen will most likely go extinct, leading to a CER that is equal (or very close) to zero even if seasonal outbreaks can occur. This renders the CER uninformative in those scenarios. We therefore devise an alternative approach for inferring outbreak risks for seasonal pathogens (involving calculating the Threshold Epidemic Risk; TER). Estimation of the TER involves calculating the probability that introduced cases will initiate a local outbreak in which a threshold number of cumulative infections is exceeded before outbreak extinction. For simple seasonal epidemic models, such as the stochastic Susceptible-Infectious-Removed model, the TER can be calculated numerically (without model simulations). For more complex models, such as stochastic host-vector models, the TER can be estimated using model simulations. We demonstrate the application of our approach by considering chikungunya virus in northern Italy as a case study. In that context, transmission is most likely in summer, when environmental conditions promote vector abundance. We show that the TER provides more useful assessments of outbreak risks than the CER, enabling practically relevant risk quantification for seasonal pathogens.

Using steady-state formula to estimate time-dependent parameters of stochastic gene transcription models

When studying stochastic gene transcription, it is important to understand how system parameters are temporally modulated in response to varying environments. Experimentally, the dynamic distribution data of RNA copy numbers measured at multiple time points are often fitted to stochastic transcription models to estimate time-dependent parameters. However, current methods require determining which parameters are time-dependent, as well as their analytical formulas, before the optimal fit. In this study, we developed a method to estimate time-dependent parameters in a classical two-state model without prior assumptions regarding the system parameters. At each measured time point, the method fitted the dynamic distribution data using a steady-state distribution formula, in which the estimated constant parameters were approximated as time-dependent parameter values at the measured time point. The accuracy of this method can be guaranteed for RNA molecules with relatively high degradation rates and genes with relatively slow responses to induction. We quantify the accuracy of the method and implemented this method on two sets of dynamic distribution data from prokaryotic and eukaryotic cells, and revealed the temporal modulation of transcription burst size in response to environmental changes.

The chemical Langevin equation for biochemical systems in dynamic environments

Modelling and simulation of complex biochemical reaction networks form cornerstones of modern biophysics. Many of the approaches developed so far capture temporal fluctuations due to the inherent stochasticity of the biophysical processes, referred to as intrinsic noise. Stochastic fluctuations, however, predominantly stem from the interplay of the network with many other - and mostly unknown - fluctuating processes, as well as with various random signals arising from the extracellular world; these sources contribute extrinsic noise. Here we provide a computational simulation method to probe the stochastic dynamics of biochemical systems subject to both intrinsic and extrinsic noise. We develop an extrinsic chemical Langevin equation-a physically motivated extension of the chemical Langevin equation- to model intrinsically noisy reaction networks embedded in a stochastically fluctuating environment. The extrinsic CLE is a continuous approximation to the Chemical Master Equation (CME) with time-varying propensities. In our approach, noise is incorporated at the level of the CME, and can account for the full dynamics of the exogenous noise process, irrespective of timescales and their mismatches. We show that our method accurately captures the first two moments of the stationary probability density when compared with exact stochastic simulation methods, while reducing the computational runtime by several orders of magnitude. Our approach provides a method that is practical, computationally efficient and physically accurate to study systems that are simultaneously subject to a variety of noise sources.

A direct comparison of methods for assessing the threat from emerging infectious diseases in seasonally varying environments

Seasonal variations in environmental conditions lead to changing infectious disease epidemic risks at different times of year. The probability that early cases initiate a major epidemic depends on the season in which the pathogen enters the population. The instantaneous epidemic risk (IER) can be tracked. This quantity is straightforward to calculate, and corresponds to the probability of a major epidemic starting from a single case introduced at time t=t₀, assuming that environmental conditions remain identical from that time onwards (i.e. for all t≥t₀). However, the threat when a pathogen enters the population in fact depends on changes in environmental conditions occurring within the timescale of the initial phase of the outbreak. For that reason, we compare the IER with a different metric: the case epidemic risk (CER). The CER corresponds to the probability of a major epidemic starting from a single case entering the population at time t=t₀, accounting for changes in environmental conditions after that time. We show how the IER and CER can be calculated using different epidemiological models (the stochastic Susceptible-Infectious-Removed model and a stochastic host-vector model that is parameterised using temperature data for Miami) in which transmission parameters vary temporally. While the IER is always easy to calculate numerically, the adaptable method we provide for calculating the CER for the host-vector model can also be applied easily and solved using widely available software tools. In line with previous research, we demonstrate that if a pathogen is likely to either invade the population or fade out on a fast timescale compared to changes in environmental conditions, the IER closely matches the CER. However, if this is not the case, the IER and the CER can be significantly different, and so the CER should be used. This demonstrates the need to consider future changes in environmental conditions carefully when assessing the risk posed by emerging pathogens.

Introduction of the Effective Photon Concentration Variable for Studying the Mechanism of Crystal Violet Photodegradation

In this work, we have proposed a new approach to study the mechanism of crystal violet (CV) photodegradation on TiO₂ surface using kinetic Monte Carlo simulation. The TiO₂ surface was considered as a set of reactive centers, which is essential in dye photodegradation. A new variable "the effective photon concentration" (I_eff ) is defined. A detailed chemical understanding of the photocatalytic reaction is provided. This approach provides a simple and effective method to find the optimal conditions of the studied system. This goal was achieved by investigating the effects of some operational parameters, including initial concentration of CV, pH, loading TiO₂ , light intensity and volume, on the degradation percent, and also, on the effective photon concentration. The perfect agreements between the experimental and simulated data at different conditions confirmed the proposed approach for describing the CV photodestruction. Also, the simulation results indicated that: (1) a significant fraction of the scattered UV irradiation into the reaction vessel does not lead to charge carrier generation; (2) the generation and recombination of charge carriers have crucial roles in the photodegradation. This is the first time that a method based on the reactive centers is employed to investigate the dye degradation by a photocatalyst.

Approximating gene transcription dynamics using steady-state formulas

Understanding how genes in a single cell respond to dynamically changing signals has been a central question in stochastic gene transcription research. Recent studies have generated massive steady-state or snapshot mRNA distribution data of individual cells, and inferred a large spectrum of kinetic transcription parameters under varying conditions. However, there have been few algorithms to convert these static data into the temporal variation of kinetic rates. Real-time imaging has been developed to monitor stochastic transcription processes at the single-cell level, but the immense technicality has prevented its application to most endogenous loci in mammalian cells. In this article, we introduced a stochastic gene transcription model with variable kinetic rates induced by unstable cellular conditions. We approximated the transcription dynamics using easily obtained steady-state formulas in the model. We tested the approximation against experimental data in both prokaryotic and eukaryotic cells and further solidified the conditions that guarantee the robustness of the method. The method can be easily implemented to provide convenient tools for quantifying dynamic kinetics and mechanisms underlying the widespread static transcription data, and may shed a light on circumventing the limitation of current bursting data on transcriptional real-time imaging.

Adapt or Perish: Evolutionary Rescue in a Gradually Deteriorating Environment

We investigate the evolutionary rescue of a microbial population in a gradually deteriorating environment, through a combination of analytical calculations and stochastic simulations. We consider a population destined for extinction in the absence of mutants, which can survive only if mutants sufficiently adapted to the new environment arise and fix. We show that mutants that appear later during the environment deterioration have a higher probability to fix. The rescue probability of the population increases with a sigmoidal shape when the product of the carrying capacity and of the mutation probability increases. Furthermore, we find that rescue becomes more likely for smaller population sizes and/or mutation probabilities if the environment degradation is slower, which illustrates the key impact of the rapidity of environment degradation on the fate of a population. We also show that our main conclusions are robust across various types of adaptive mutants, including specialist and generalist ones, as well as mutants modeling antimicrobial resistance evolution. We further express the average time of appearance of the mutants that do rescue the population and the average extinction time of those that do not. Our methods can be applied to other situations with continuously variable fitnesses and population sizes, and our analytical predictions are valid in the weak-to-moderate mutation regime.

Is mammography screening beneficial: An individual-based stochastic model for breast cancer incidence and mortality

The benefits of mammography screening have been controversial, with conflicting findings from various studies. We hypothesize that unmeasured heterogeneity in tumor aggressiveness underlies these conflicting results. Based on published data from the Canadian National Breast Screening Study (CNBSS), we develop and parameterize an individual-based mechanistic model for breast cancer incidence and mortality that tracks five stages of breast cancer progression and incorporates the effects of age on breast cancer incidence and all-cause mortality. The model accurately reproduces the reported outcomes of the CNBSS. By varying parameters, we predict that the benefits of mammography depend on the effectiveness of cancer treatment and tumor aggressiveness. In particular, patients with the most rapidly growing or potentially largest tumors have the highest benefit and least harm from the screening, with only a relatively small effect of age. However, the model predicts that confining mammography to populations with a high risk of acquiring breast cancer increases the screening benefit only slightly compared with the full population.

Stochastic simulation algorithms for computational systems biology: Exact, approximate, and hybrid methods

Nowadays, mathematical modeling is playing a key role in many different research fields. In the context of system biology, mathematical models and their associated computer simulations constitute essential tools of investigation. Among the others, they provide a way to systematically analyze systems perturbations, develop hypotheses to guide the design of new experimental tests, and ultimately assess the suitability of specific molecules as novel therapeutic targets. To these purposes, stochastic simulation algorithms (SSAs) have been introduced for numerically simulating the time evolution of a well-stirred chemically reacting system by taking proper account of the randomness inherent in such a system. In this work, we review the main SSAs that have been introduced in the context of exact, approximate, and hybrid stochastic simulation. Specifically, we will introduce the direct method (DM), the first reaction method (FRM), the next reaction method (NRM) and the rejection-based SSA (RSSA) in the area of exact stochastic simulation. We will then present the τ-leaping method and the chemical Langevin method in the area of approximate stochastic simulation and an implementation of the hybrid RSSA (HRSSA) in the context of hybrid stochastic-deterministic simulation. Finally, we will consider the model of the sphingolipid metabolism to provide an example of application of SSA to computational system biology by exemplifying how different simulation strategies may unveil different insights into the investigated biological phenomenon. This article is categorized under: Models of Systems Properties and Processes > Mechanistic Models Analytical and Computational Methods > Computational Methods.

Efficient anticorrelated variance reduction for stochastic simulation of biochemical reactions

We investigate the computational challenge of improving the accuracy of the stochastic simulation estimation by inducing negative correlation through the anticorrelated variance reduction technique. A direct application of the technique to the stochastic simulation algorithm (SSA), employing the inverse transformation, is not efficient for simulating large networks because its computational cost is similar to the sum of independent simulation runs. We propose in this study a new algorithm that employs the propensity bounds of reactions, introduced recently in their rejection-based SSA, to correlate and synchronise the trajectories during the simulation. The selection of reaction firings by our approach is exact due to the rejection-based mechanism. In addition, by applying the anticorrelated variance technique to select reaction firings, our approach can induce substantial correlation between realisations, hence reducing the variance of the estimator. The computational advantage of our rejection-based approach in comparison with the traditional inverse transformation is that it only needs to maintain a single data structure storing propensity bounds of reactions, which is updated infrequently, hence achieving better performance.

Efficient finite-difference method for computing sensitivities of biochemical reactions

Sensitivity analysis of biochemical reactions aims at quantifying the dependence of the reaction dynamics on the reaction rates. The computation of the parameter sensitivities, however, poses many computational challenges when taking stochastic noise into account. This paper proposes a new finite-difference method for efficiently computing sensitivities of biochemical reactions. We employ propensity bounds of reactions to couple the simulation of the nominal and perturbed processes. The exactness of the simulation is preserved by applying the rejection-based mechanism. For each simulation step, the nominal and perturbed processes under our coupling strategy are synchronized and often jump together, increasing their positive correlation and hence reducing the variance of the estimator. The distinctive feature of our approach in comparison with existing coupling approaches is that it only needs to maintain a single data structure storing propensity bounds of reactions during the simulation of the nominal and perturbed processes. Our approach allows to compute sensitivities of many reaction rates simultaneously. Moreover, the data structure does not require to be updated frequently, hence improving the computational cost. This feature is especially useful when applied to large reaction networks. We benchmark our method on biological reaction models to prove its applicability and efficiency.

Incorporating extrinsic noise into the stochastic simulation of biochemical reactions: A comparison of approaches

The stochastic simulation algorithm (SSA) has been widely used for simulating biochemical reaction networks. SSA is able to capture the inherently intrinsic noise of the biological system, which is due to the discreteness of species population and to the randomness of their reciprocal interactions. However, SSA does not consider other sources of heterogeneity in biochemical reaction systems, which are referred to as extrinsic noise. Here, we extend two simulation approaches, namely, the integration-based method and the rejection-based method, to take extrinsic noise into account by allowing the reaction propensities to vary in time and state dependent manner. For both methods, new efficient implementations are introduced and their efficiency and applicability to biological models are investigated. Our numerical results suggest that the rejection-based method performs better than the integration-based method when the extrinsic noise is considered.

A Critical Comparison of Rejection-Based Algorithms for Simulation of Large Biochemical Reaction Networks

The rejection-based simulation technique has been applying to improve the computational efficiency of the stochastic simulation algorithm (SSA) in simulating large reaction networks, which are required for a thorough understanding of biological systems. We compare two recently proposed simulation methods, namely the composition-rejection algorithm (SSA-CR) and the rejection-based SSA (RSSA), aiming for this purpose. We discuss the right interpretation of the rejection-based technique used in these algorithms in order to make an informed choice when dealing with different aspects of biochemical networks. We provide the theoretical analysis as well as the detailed runtime comparison of these algorithms on concrete biological models. We highlight important factors that are omitted in previous analysis of these algorithms. The numerical comparison shows that for reaction networks where the search cost is expensive then SSA-CR is more efficient, and for reaction networks where the update cost is dominant, often the case in practice, then RSSA should be the choice.

Stochastic simulation of biochemical reactions with partial-propensity and rejection-based approaches

We present in this paper a new exact algorithm for improving performance of exact stochastic simulation algorithm. The algorithm is developed on concepts of the partial-propensity and the rejection-based approaches. It factorizes the propensity bounds of reactions and groups factors by common reactant species for selecting next reaction firings. Our algorithm provides favorable computational advantages for simulating of biochemical reaction networks by reducing the cost for selecting the next reaction firing to scale with the number of chemical species and avoiding expensive propensity updates during the simulation. We present the details of our new algorithm and benchmark it on concrete biological models to demonstrate its applicability and efficiency.

Efficient Constant-Time Complexity Algorithm for Stochastic Simulation of Large Reaction Networks

Exact stochastic simulation is an indispensable tool for a quantitative study of biochemical reaction networks. The simulation realizes the time evolution of the model by randomly choosing a reaction to fire and update the system state according to a probability that is proportional to the reaction propensity. Two computationally expensive tasks in simulating large biochemical networks are the selection of next reaction firings and the update of reaction propensities due to state changes. We present in this work a new exact algorithm to optimize both of these simulation bottlenecks. Our algorithm employs the composition-rejection on the propensity bounds of reactions to select the next reaction firing. The selection of next reaction firings is independent of the number reactions while the update of propensities is skipped and performed only when necessary. It therefore provides a favorable scaling for the computational complexity in simulating large reaction networks. We benchmark our new algorithm with the state of the art algorithms available in literature to demonstrate its applicability and efficiency.

A mechanistic stochastic framework for regulating bacterial cell division

How exponentially growing cells maintain size homeostasis is an important fundamental problem. Recent single-cell studies in prokaryotes have uncovered the adder principle, where cells add a fixed size (volume) from birth to division, irrespective of their size at birth. To mechanistically explain the adder principle, we consider a timekeeper protein that begins to get stochastically expressed after cell birth at a rate proportional to the volume. Cell-division time is formulated as the first-passage time for protein copy numbers to hit a fixed threshold. Consistent with data, the model predicts that the noise in division timing increases with size at birth. Intriguingly, our results show that the distribution of the volume added between successive cell-division events is independent of the newborn cell size. This was dramatically seen in experimental studies, where histograms of the added volume corresponding to different newborn sizes collapsed on top of each other. The model provides further insights consistent with experimental observations: the distribution of the added volume when scaled by its mean becomes invariant of the growth rate. In summary, our simple yet elegant model explains key experimental findings and suggests a mechanism for regulating both the mean and fluctuations in cell-division timing for controlling size.

Stochastic Simulation of Biomolecular Networks in Dynamic Environments

Simulation of biomolecular networks is now indispensable for studying biological systems, from small reaction networks to large ensembles of cells. Here we present a novel approach for stochastic simulation of networks embedded in the dynamic environment of the cell and its surroundings. We thus sample trajectories of the stochastic process described by the chemical master equation with time-varying propensities. A comparative analysis shows that existing approaches can either fail dramatically, or else can impose impractical computational burdens due to numerical integration of reaction propensities, especially when cell ensembles are studied. Here we introduce the Extrande method which, given a simulated time course of dynamic network inputs, provides a conditionally exact and several orders-of-magnitude faster simulation solution. The new approach makes it feasible to demonstrate—using decision-making by a large population of quorum sensing bacteria—that robustness to fluctuations from upstream signaling places strong constraints on the design of networks determining cell fate. Our approach has the potential to significantly advance both understanding of molecular systems biology and design of synthetic circuits.

Simulation algorithms have become indispensable tools in modern quantitative biology, providing deep insight into many biochemical systems, including gene regulatory networks. However, current stochastic simulation approaches handle the effects of fluctuating extracellular signals and upstream processes poorly, either failing to give qualitatively reliable predictions or being very inefficient computationally. Here we introduce the Extrande method, a novel approach for simulation of biomolecular networks embedded in the dynamic environment of the cell and its surroundings. The method is accurate and computationally efficient, and hence fills an important gap in the field of stochastic simulation. In particular, we employ it to study a bacterial decision-making network and demonstrate that robustness to fluctuations from upstream signaling places strong constraints on the design of networks determining cell fate.