Stochastic Simulation Algorithm Research Articles

Synchronization and Random Attractors in Reaction Jump Processes

AbstractThis work explores a synchronization-like phenomenon induced by common noise for continuous-time Markov jump processes given by chemical reaction networks. Based on Gillespie’s stochastic simulation algorithm, a corresponding random dynamical system is formulated in a two-step procedure, at first for the states of the embedded discrete-time Markov chain and then for the augmented Markov chain including random jump times. We uncover a time-shifted synchronization in the sense that—after some initial waiting time—one trajectory exactly replicates another one with a certain time delay. Whether or not such a synchronization behavior occurs depends on the combination of the initial states. We prove this partial time-shifted synchronization for the special setting of a birth-death process by analyzing the corresponding two-point motion of the embedded Markov chain and determine the structure of the associated random attractor. In this context, we also provide general results on existence and form of random attractors for discrete-time, discrete-space random dynamical systems.

Time-dependent earthquake-fire coupling fragility analysis under limited prior knowledge: A perspective from type-2 fuzzy probability

Earthquake-triggered fire domino scenarios (E-FDSs) arise frequently from the interaction between earthquakes and chemical installations, resulting in catastrophic multi-hazard coupling events. The complicated mutually amplified phenomena between natural disasters and chemical accidents significantly aggravates the escalation of domino accidents, which has posed great challenges for modeling and preventing E-FDSs. Under this impetus, this work proposes an advanced type-2 fuzzy probabilistic methodology to obtain the time-dependent failure probability of steel cylindrical tanks (SCTs) subjected to the earthquake-fire sequence. To cope with the limited prior knowledge on E-FDSs, a basic universal is established to describe the fire resistance attenuation caused by the seismic damage. The coupling failure criterion of SCTs is formulated by a type-2 fuzzy time-dependent limit state equation. A credibility-based stochastic simulation algorithm is developed for the hybrid uncertainty analysis (combining ambiguity and stochasticity). The proposed methodology is validated by case studies of a 5000 m3 fixed roof tank. Compared to the existing accident probability model, the proposed methodology can not only capture the fire resistance attenuation caused by the seismic damage but also provide a dynamic estimation of tank failure probability with respect to the fire exposure time. The proposed methodology can effectively and dynamically capture the accident evolution process, which in turn helps mitigate and prevent the spatiotemporal propagation of domino effects.

Genetic Network Design Automation with LOICA.

Genetic design automation (GDA) is the use of computer-aided design (CAD) in designing genetic networks. GDA tools are necessary to create more complex synthetic genetic networks in a high-throughput fashion. At the core of these tools is the abstraction of a hierarchy of standardized components. The components' input, output, and interactions must be captured and parametrized from relevant experimental data. Simulations of genetic networks should use those parameters and include the experimental context to be compared with the experimental results.This chapter introduces Logical Operators for Integrated Cell Algorithms (LOICA), a Python package used for designing, modeling, and characterizing genetic networks using a simple object-oriented design abstraction. LOICA represents different biological and experimental components as classes that interact to generate models. These models can be parametrized by direct connection to the Flapjack experimental data management platform to characterize abstracted components with experimental data. The models can be simulated using stochastic simulation algorithms or ordinary differential equations with varying noise levels. The simulated data can be managed and published using Flapjack alongside experimental data for comparison. LOICA genetic network designs can be represented as graphs and plotted as networks for visual inspection and serialized as Python objects or in the Synthetic Biology Open Language (SBOL) format for sharing and use in other designs.

A numerical methodology for estimating site-specific cascading earthquake and tsunami dynamic loading on critical infrastructure

A numerical methodology for deriving cascading earthquake and tsunami loading patterns for use in critical infrastructure design is proposed. Earthquake source parameters are first used to generate synthetic ground motion acceleration time series at the site of a critical infrastructure using an open-source stochastic extended finite-source simulation algorithm. Using the same earthquake source information, the tsunami phases of generation, propagation, and inundation are modeled using a coupling approach for determining tsunami wave forces. The coupling approach combines (1) a code based on non-linear shallow water equations solved by the Finite Volume method, which is used to capture the tsunami flow characteristics from generation through the inundation phase until relatively close to the structure, coupled with (2) an open-source computational fluid simulation code, which relies on solving the Navier–Stokes equations using the Smoothed Particle Hydrodynamics method. An appendix presents the validation of the numerical methodology through correlation of the numerical simulations with field data recorded during the 2011 Japan Earthquake and Tsunami. The methodology is then employed to characterize the ground shaking and tsunami loading patterns at the site of the container terminal of the Sines deep-water seaport in Portugal, when subjected to the historic 1755 Great Lisbon Earthquake and Tsunami plausible scenarios. For the site, the generated ground motion acceleration time series show peak accelerations exceeding Eurocode 8 design values. Tsunami flow depths and velocities were computed in horizontal and vertical directions and could be used in the assessment of the existing infrastructure or the design of the future expansion of the port.

Stochastic Simulation Algorithms for Iterative Solution of the Lamé Equation

Stochastic Simulation Algorithms for Iterative Solution of the Lamé Equation

IDESS: a toolbox for identification and automated design of stochastic gene circuits.

One of the main causes hampering predictability during the model identification and automated design of gene circuits in synthetic biology is the effect of molecular noise. Stochasticity may significantly impact the dynamics and function of gene circuits, specially in bacteria and yeast due to low mRNA copy numbers. Standard stochastic simulation methods are too computationally costly in realistic scenarios to be applied to optimization-based design or parameter estimation. In this work, we present IDESS (Identification and automated Design of Stochastic gene circuitS), a software toolbox for optimization-based design and model identification of gene regulatory circuits in the stochastic regime. This software incorporates an efficient approximation of the Chemical Master Equation as well as a stochastic simulation algorithm-both with GPU and CPU implementations-combined with global optimization algorithms capable of solving Mixed Integer Nonlinear Programming problems. The toolbox efficiently addresses two types of problems relevant in systems and synthetic biology: the automated design of stochastic synthetic gene circuits, and the parameter estimation for model identification of stochastic gene regulatory networks. IDESS runs under the MATLAB environment and it is available under GPLv3 license at https://doi.org/10.5281/zenodo.7788692.

Stochastic simulation of electron transport in a strong electrical field in low-dimensional heterostructures

Abstract In this paper we develop a stochastic simulation algorithm for electron transport in a DA-pHEMT heterostructure. Mathematical formulation of the problem of electron gas transport in the heterostructure in the form of a coupled system of Poisson, Schrödinger and kinetic Boltzmann equations is given. A Monte Carlo model of electron transport in DA-pHEMT heterostructures which accounts for multivalley parabolic band structure, as well as relevant formulas for calculating electron scattering rates and scattering phase functions on polar optical, intervalley phonons and on impurities are developed. The results of a computational experiment involving the solution of the system of Poisson–Schrödinger–Boltzmann equations for the AlGaAs/GaAs/InGaAs/GaAs/AlGaAs heterostructure are presented. The distribution of electrons by energy subband in the main and satellite valleys and the field dependences of the electron drift velocity in each valley are calculated. It was discovered that there is no spatial transfer of electrons into wide-gap AlGaAs layers due to high barriers created by modulated-doped impurities. A comparative analysis of the electron drift velocities in the studied DA-pHEMT heterostructures and in the unstrained layer of the InGaAs is given.

From Compositional Petri Net Modeling to Macro and Micro Simulation by Means of Stochastic Simulation and Agent-Based Models

Computational modeling has become a widespread approach for studying real-world phenomena by using different modeling perspectives, in particular, the microscopic point of view concentrates on the behavior of the single components and their interactions from which the global system evolution emerges, while the macroscopic point of view represents the system’s overall behavior abstracting as much as possible from that of the single components. The preferred point of view depends on the effort required to develop the model, on the detail level of the available information about the system to be modeled, and on the type of measures that are of interest to the modeler; each point of view may lead to a different modeling language and simulation paradigm. An approach adequate for the microscopic point of view is Agent-Based Modeling and Simulation, which has gained popularity in the last few decades but lacks a formal definition common to the different tools supporting it. This may lead to modeling mistakes and wrong interpretation of the results, especially when comparing models of the same system developed according to different points of view. The aim of the work described in this paper is to provide a common compositional modeling language from which both a macro and a micro simulation model can be automatically derived: these models are coherent by construction and may be studied through different simulation approaches and tools. A framework is thus proposed in which a model can be composed using a Petri Net formalism and then studied through both an Agent-Based Simulation and a classical Stochastic Simulation Algorithm, depending on the study goal.

Optimal explicit stabilized postprocessed τ-leap method for the simulation of chemical kinetics

The simulation of chemical kinetics involving multiple scales constitutes a modeling challenge (from ordinary differential equations to Markov chains) and a computational challenge (multiple scales, large dynamical systems, time step restrictions). In this paper, we propose a new discrete stochastic simulation algorithm: the postprocessed second kind stabilized orthogonal τ-leap Runge–Kutta method (PSK-τ-ROCK). In the context of chemical kinetics, this method can be seen as a stabilization of Gillespie's explicit τ-leap algorithm combined with a postprocessor. The stabilized procedure allows to simulate problems with multiple scales (stiff), while the postprocessing procedure allows to approximate the invariant measure (e.g. mean and variance) of ergodic stochastic dynamical systems. We prove stability and accuracy of the PSK-τ-ROCK for a reference system. Numerical experiments illustrate the high reliability and efficiency of the scheme when compared to other τ-leap methods.

The regime-conversion method: a hybrid technique for simulating well-mixed chemical reaction networks

There exist several methods for simulating biological and physical systems as represented by chemical reaction networks. Systems with low numbers of particles are frequently modeled as discrete-state Markov jump processes and are typically simulated via a stochastic simulation algorithm (SSA). An SSA, while accurate, is often unsuitable for systems with large numbers of individuals, and can become prohibitively expensive with increasing reaction frequency. Large systems are often modeled deterministically using ordinary differential equations, sacrificing accuracy and stochasticity for computational efficiency and analytical tractability. In this paper, we present a novel hybrid technique for the accurate and efficient simulation of large chemical reaction networks. This technique, which we name the regime-conversion method, couples a discrete-state Markov jump process to a system of ordinary differential equations by simulating a reaction network using both techniques simultaneously. Individual molecules in the network are represented by exactly one regime at any given time, and may switch their governing regime depending on particle density. In this manner, we model high copy-number species using the cheaper continuum method and low copy-number species using the more expensive, discrete-state stochastic method to preserve the impact of stochastic fluctuations at low copy number. The motivation, as with similar methods, is to retain the advantages while mitigating the shortfalls of each method. We demonstrate the performance and accuracy of our method for several test problems that exhibit varying degrees of inter-connectivity and complexity by comparing averaged trajectories obtained from both our method and from exact stochastic simulation.

An Implementation of The Method of Moments on Chemical Systems with Constant and Time-dependent Rates

Among numerical techniques used to facilitate the analysis of biochemical reactions, we can use the method of moments to directly approximate statistics such as the mean numbers of molecules. The method is computationally viable in time and memory, compared to solving the chemical master equation (CME) which is notoriously expensive. In this study, we apply the method of moments to a chemical system with a constant rate representing a vascular endothelial growth factor (VEGF) model, as well as another system with time-dependent propensities representing the susceptible, infected, and recovered (SIR) model with periodic contact rate. We assess the accuracy of the method using comparisons with approximations obtained by the stochastic simulation algorithm (SSA) and the chemical Langevin equation (CLE). The VEGF model is of interest because of the role of VEGF in the growth of cancer and other inflammatory diseases and the potential use of anti-VEGF therapies in the treatment of cancer. The SIR model is a popular epidemiological model used in studying the spread of various infectious diseases in a population.

Tensor product approach to modelling epidemics on networks

To improve mathematical models of epidemics it is essential to move beyond the traditional assumption of homogeneous well–mixed population and involve more precise information on the network of contacts and transport links by which a stochastic process of the epidemics spreads. In general, the number of states of the network grows exponentially with its size, and a master equation description suffers from the curse of dimensionality. Almost all methods widely used in practice are versions of the stochastic simulation algorithm (SSA), which is notoriously known for its slow convergence. In this paper we numerically solve the chemical master equation for an SIR model on a general network using recently proposed tensor product algorithms. In numerical experiments we show that tensor product algorithms converge much faster than SSA and deliver more accurate results, which becomes particularly important for uncovering the probabilities of rare events, e.g. for number of infected people to exceed a (high) threshold.

Retracted: An Evaluation Approach for Physical Education Teaching Practice Quality Using Stochastic Simulation Algorithm

Retracted: An Evaluation Approach for Physical Education Teaching Practice Quality Using Stochastic Simulation Algorithm

Retracted: A Study on the Usefulness of Stochastic Simulation Algorithms for Teaching and Learning in College Physical Education Classrooms

Retracted: A Study on the Usefulness of Stochastic Simulation Algorithms for Teaching and Learning in College Physical Education Classrooms

Disruption in the regulation of casein kinase 2 in circadian rhythm leads to pathological states: cancer, diabetes and neurodegenerative disorders.

Circadian rhythm maintains the sleep-wake cycle in biological systems. Various biological activities are regulated and modulated by the circadian rhythm, disruption of which can result in onset of diseases. Robust rhythms of phosphorylation profiles and abundances of PERIOD (PER) proteins are thought to be the master keys that drive circadian clock functions. The role of casein kinase 2 (CK2) in circadian rhythm via its direct interactions with the PER protein has been extensively studied; however, the exact mechanism by which it affects circadian rhythms at the molecular level is not known. Here, we propose an extended circadian rhythm model in Drosophila that incorporates the crosstalk between the PER protein and CK2. We studied the regulatory role of CK2 in the dynamics of PER proteins involved in circadian rhythm using the stochastic simulation algorithm. We observed that variations in the concentration of CK2 in the circadian rhythm model modulates the PER protein dynamics at different cellular states, namely, active, weakly active, and rhythmic death. These oscillatory states may correspond to distinct pathological cellular states of the living system. We find molecular noise at the expression level of CK2 to switch normal circadian rhythm to any of the three above-mentioned circadian oscillatory states. Our results suggest that the concentration levels of CK2 in the system has a strong impact on its dynamics, which is reflected in the time evolution of PER protein. We believe that our findings can contribute towards understanding the molecular mechanisms of circadian dysregulation in pathways driven by the PER mutant genes and their pathological states, including cancer, obesity, diabetes, neurodegenerative disorders, and socio-psychological disease.

The impossible challenge of estimating non-existent moments of the Chemical Master Equation.

The Chemical Master Equation (CME) is a set of linear differential equations that describes the evolution of the probability distribution on all possible configurations of a (bio-)chemical reaction system. Since the number of configurations and therefore the dimension of the CME rapidly increases with the number of molecules, its applicability is restricted to small systems. A widely applied remedy for this challenge is moment-based approaches which consider the evolution of the first few moments of the distribution as summary statistics for the complete distribution. Here, we investigate the performance of two moment-estimation methods for reaction systems whose equilibrium distributions encounter fat-tailedness and do not possess statistical moments. We show that estimation via stochastic simulation algorithm (SSA) trajectories lose consistency over time and estimated moment values span a wide range of values even for large sample sizes. In comparison, the method of moments returns smooth moment estimates but is not able to indicate the non-existence of the allegedly predicted moments. We furthermore analyze the negative effect of a CME solution's fat-tailedness on SSA run times and explain inherent difficulties. While moment-estimation techniques are a commonly applied tool in the simulation of (bio-)chemical reaction networks, we conclude that they should be used with care, as neither the system definition nor the moment-estimation techniques themselves reliably indicate the potential fat-tailedness of the CME's solution.

Assessment of the Spatial Variability and Uncertainty of Shreddable Pruning Biomass in an Olive Grove Based on Canopy Volume and Tree Projected Area

Olive pruning residues are a by-product that can be applied to soil or used for energy production in a circular economy model. Its benefits depend on the amount of pruning, which varies greatly within farms. This study aimed to investigate the spatial variability of shreddable olive pruning in a traditional olive grove in Córdoba (Spain) with an area of 15 ha and trees distanced 12.5 m from each other. To model the spatial variability of shreddable olive pruning, geostatistical methods of stochastic simulation were applied to three correlated variables measured on sampled trees: the crown projected area (n = 928 trees), the crown volume (n = 167) and the amount of shreddable pruning (n = 59). Pearson’s correlation between pairs of variables varied from 0.71 to 0.76. The amount of pruning showed great variability, ranging from 7.6 to 76 kg tree−1, with a mean value of 37 kg tree−1. Using exponential and spherical variogram models, the spatial continuity of the variables under study was established. Shreddable dry pruning weight values showed spatial autocorrelation up to 180 m. The spatial uncertainty of the estimation was obtained using sequential simulation algorithms. Stochastic simulation algorithms provided 150 possible images of the amount of shreddable pruning on the farm, using tree projected area and crown volume as secondary information. The interquartile range and 90% prediction interval were used as indicators of the uncertainty around the mean value. Uncertainty validation was performed using accuracy plots and the associated G-statistic. Results indicate with high confidence (i.e., low uncertainty) that shreddable dry pruning weight in the mid-western area of the farm will be much lower than the rest of the farm. In the same way, results show with high confidence that dry pruning weight will be much higher in a small area in the middle east of the farm. The values of the G-statistic ranged between 0.89 and 0.90 in the tests performed. The joint use of crown volume and projected areas is valuable in estimating the spatial variability of the amount of pruning. The study shows that the use of prediction intervals enables the evaluation of farm areas and informed management decisions with a low level of risk. The methodology proposed in this work can be extrapolated to other 3D crops without requiring modifications. On a larger scale, it can be useful for predicting optimal locations for biomass plants, areas with high potential as carbon sinks or areas requiring special soil protection measures.

Optimization of declared capacity for electrical fused magnesia furnaces participating in primary frequency control

In the frequency response of the new power system, fused magnesium as a load side resource has great potential for frequency modulation. Due to the need to adapt to the medium and long term power market demand, the existing fused magnesium simulation methods which only take single ton energy consumption as the index are difficult to determine the declared capacity. Therefore, according to the requirement of annual time scale and the randomness of molten magnesium production process, a stochastic production simulation algorithm framework was designed to simulate the operation state of molten magnesium load. Based on the Markov chain Monte Carlo method and the idea of integer random distribution, the transition and duration of each state were simulated. According to the rule of two-step system, an optimization model for the declared capacity of fused magnesium was established, and the optimal declared capacity and enterprise income were solved by the ergodic optimization algorithm, so as to achieve the goal of frequency response on the load side. The simulation example shows that the proposed method can well characterize the production process of electrical fused magnesium load, calculate the optimal declared capacity for the scenario of power loss in the medium and long term market, which brings economic benefits to enterprises and provides fast frequency response resources for power grid.

Two stochastic algorithms for solving elastostatics problems governed by the Lamé equation

Abstract In this paper, we construct stochastic simulation algorithms for solving an elastostatics problem governed by the Lamé equation. Two different stochastic simulation methods are suggested: (1) a method based on a random walk on spheres, which is iteratively applied to anisotropic diffusion equations that are related through the mixed second-order derivatives (this method is meshless and can be applied to boundary value problems for complicated domains); (2) a randomized algorithm for solving large systems of linear algebraic equations that is the core of this method. It needs a mesh formation, but even for very fine grids, the algorithm shows a high efficiency. Both methods are scalable and can be easily parallelized.

Stochastic Simulations of Casual Groups

Free-forming or casual groups are groups in which individuals are in face-to-face interactions and are free to maintain or terminate contact with one another, such as clusters of people at a cocktail party, play groups in a children’s playground or shopping groups in a mall. Stochastic models of casual groups assume that group sizes are the products of natural processes by which groups acquire and lose members. The size distributions predicted by these models have been the object of controversy since their derivation in the 1960s because of the neglect of fluctuations around the mean values of random variables that characterize a collection of groups. Here, we check the validity of these mean-field approximations using an exact stochastic simulation algorithm to study the processes of the acquisition and loss of group members. In addition, we consider the situation where the appeal of a group of size i to isolates is proportional to iα. We find that, for α≤1, the mean-field approximation fits the equilibrium simulation results very well, even for a relatively small population size N. However, for α>1, this approximation scheme fails to provide a coherent description of the distribution of group sizes. We find a discontinuous phase transition at αc>1 that separates the regime where the variance of the group size does not depend on N from the regime where it grows linearly with N. In the latter regime, the system is composed of a single large group that coexists with a large number of isolates. Hence, the same underlying acquisition-and-loss process can explain the existence of small, temporary casual groups and of large, stable social groups.