Stochastic Simulation Algorithm Research Articles

Development and Application of a Tornado Database for the Chinese Mainland

This study develops a preliminary tornado database for the Chinese mainland using information provided by the Yearbook of Meteorological Disasters in China (YMDC) from 2003 to 2019 and other public media report data. This database includes tornado occurrence time, geographical location, intensity, and damage descriptions. A modified Enhanced Fujita (EF) scale criterion adapted to the description of the damage indicator (DI) and degree of damage (DOD) in the YMDC and media reports is presented to better estimate the tornado intensity. The spatial and temporal distribution characteristics of tornadoes in Chinese mainland are examined. A stochastic simulation algorithm is proposed to perform a risk assessment of tornado hazards. The occurrence of tornadoes is randomly sampled using negative binomial distribution. The kernel density estimation (KDE) method based on the Gaussian kernel is applied to estimate the probability density of geographically dependent tornado occurrence before randomly generating tornado occurrence locations. The probability and conditional probability of tornado occurrence with different intensity levels are calculated for different counties in Jiangsu and Guangdong Provinces using Bayes’ theorem. The database provides a forward step toward rational assessment of tornado-induced disasters in Chinese mainland. The database can be accessed at https://doi.org/10.4121/19801633.v1 (Zhang et al., 2022).

Intraspecific predator interference promotes biodiversity in ecosystems.

Explaining biodiversity is a fundamental issue in ecology. A long-standing puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individual-based modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.

Intraspecific predator interference promotes biodiversity in ecosystems

Explaining biodiversity is a fundamental issue in ecology. A long-standing puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individual-based modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.

Self-organized biodiversity in biotic resource systems through pairwise encounters

A central issue in ecology is understanding what determines biodiversity in nature, particularly in biotic resource systems that typically lack cross-feeding. In this study, we demonstrate that by incorporating pairwise encounters among individuals of the same consumer species, a diverse array of consumer species can self-organize to coexist in a well-mixed system with one or a few biotic resource species. The coexistence modes may manifest as either stable steady states or self-organized oscillations. Notably, all coexistence states are robust to stochasticity, as confirmed by both stochastic simulation algorithms and individual-based modeling. Our model quantitatively illustrates species distribution patterns in various ecological communities and can be broadly used to explain biodiversity in many biotic resource systems.

A novel approach for fractional pharmacokinetics modeling and integrating stochastic simulation techniques using Sibuya distribution

Fractional compartment models are widely used to simulate pharmacokinetics with power-law behavior. In this work, we study an approach, previously mainly used in Epidemiology for the age model, employing the Sibuya distribution, a discrete-time, non-Markov probability distribution, to effectively emulate Mittag-Leffler functions. We demonstrate the validity of our method by applying it to various pharmacokinetic datasets and comparing the results with traditional methods like the Numerical Inverse Laplace Transform (NILT) algorithm. Additionally, we propose two stochastic simulation algorithms, where the Monte Carlo SSA allows us to capture the randomness of data and the Gillespie-like stochastic numerical simulation technique can optimize computational resources in fractional analysis.

In enzymatic reactions, the reverse reaction reduces product noise

Enzymatic reactions are essential for most cellular reactions and ubiquitous in living organisms. In the present study, we explore the pivotal role of the reverse reaction in enzymatic reactions. It is a powerful noise-buffering motif. By SSA (stochastic simulation algorithm), a remarkable 32% reduction of product CV (coefficient of variation) was observed. To better understand the causes, we split the upstream noise. The product CV reduction is more than 35% for the noise inherited from the enzyme but merely 6%–21% for that from the substrate. It implies that the system applies different strategies to different upstream noises. We identified two leading causes responsible for noise attenuation. A cell is well designed to control its intracellular noise, and to acquire wisdom from nature is always enjoyable.

Tensor product algorithms for inference of contact network from epidemiological data

We consider a problem of inferring contact network from nodal states observed during an epidemiological process. In a black-box Bayesian optimisation framework this problem reduces to a discrete likelihood optimisation over the set of possible networks. The cardinality of this set grows combinatorially with the number of network nodes, which makes this optimisation computationally challenging. For each network, its likelihood is the probability for the observed data to appear during the evolution of the epidemiological process on this network. This probability can be very small, particularly if the network is significantly different from the ground truth network, from which the observed data actually appear. A commonly used stochastic simulation algorithm struggles to recover rare events and hence to estimate small probabilities and likelihoods. In this paper we replace the stochastic simulation with solving the chemical master equation for the probabilities of all network states. Since this equation also suffers from the curse of dimensionality, we apply tensor train approximations to overcome it and enable fast and accurate computations. Numerical simulations demonstrate efficient black-box Bayesian inference of the network.

Stochastic simulation of exciton transport in semiconductor heterostructures

Abstract Stochastic simulation algorithm for solving exciton transport in a 3D layered semiconductor heterostructure is developed. The problem is governed by a transient drift-diffusion-recombination equation with Dirichlet and Neumann mixed boundary conditions. The semiconductor is represented as an infinite multilayer of finite thickness along the transverse coordinate z. The multilayer is composed by a set of sublayers of different materials so that the excitons have different diffusion and recombination coefficients in each layer. Continuity of solutions and fluxes at the plane interfaces between layers are imposed. The stochastic simulation algorithm solves the transport problem by tracking exciton trajectories in accordance with the probability distributions represented through the Green function of the problem in each sublayer. The method is meshless, the excitons jump only over the plane boundaries of the layers. This explains the high efficiency of the method. Simulation results for transport problems with different mixed boundary conditions are presented.

Stochastic Simulation Algorithm for Solving the System of Lame Equations for Two- and Three-Dimensional Domains by Combining the Slobodianskii Representation, the Method of Fundamental Solutions and a Stochastic Projection Method

Stochastic Simulation Algorithm for Solving the System of Lame Equations for Two- and Three-Dimensional Domains by Combining the Slobodianskii Representation, the Method of Fundamental Solutions and a Stochastic Projection Method

Reservoir characteristics analysis and favorable area prediction of Zhuhai Formation, Wenchang A oilfield, South China Sea

Wenchang A oilfield is a recently discovered low-permeability oilfield in the western South China Sea. The exploration target is the Zhuhai Formation, which reservoir lithology changes significantly, and the distribution of favorable reservoirs is unclear. In this study, the reservoir characteristics of the Zhuhai Formation in the Wenchang A oilfield, South China Sea, were analyzed through core identification, thin section identification, physical property testing, pre-stack geostatistical inversion, and frequency-based AVO inversion. Pre-stack geostatistical inversion is based on geostatistics and combines the seismic inversion algorithm with the stochastic sequential simulation algorithm. While the frequency-varying AVO inversion method is an extension of the conventional AVO inversion method, with its core technology being spectrum decomposition. After frequency division, multiple data sets directly participate in the calculation, which enhances the stability and accuracy of inversion.We also predict the favorable reservoir areas of six small layers in the Zhuhai Formation. The results show that 1) the reservoir lithology of the Zhuhai Formation mainly consists of feldspar quartz sandstone and feldspar quartz sandstone, indicating high compositional maturity. The porosity of the reservoir ranges from 14.0% to 19.0%, with an average value of 16.3%. The permeability of the reservoir ranges from 3.1 mD to 126.1 mD, with an average value of 22.4 mD, indicating a medium porosity and low permeability reservoir. The supporting structure of a sandstone reservoir is particle-supported, with the main contact being the “point-line” contact between particles. The main types of pores in sandstone are primary intergranular pores and secondary intergranular dissolved pores. Some pores are feldspar-dissolved pores, and occasionally there are hetero-based micropore. 2) The diagenesis of the Zhuhai Formation reservoir mainly includes compaction, cementation, and dissolution. Cementation and dissolution have minimal impact on the physical properties of reservoirs. The porosity loss of the reservoir after compaction ranges between 16.0% and 27.9%, and the compaction rate ranges between 40.0% and 69.8%, indicating a medium compaction diagenetic facies. 3) The favorable areas of each layer of the Zhuhai Formation are mainly concentrated in the southwest of the study area. The upper layer of ZH1I, the upper layer of ZH1II and the lower layer of ZH1II exhibit the best physical properties, the thickest favorable sand body, the strongest oil and gas display, the widest range of favorable areas, and the greatest exploration potential. The favorable exploration potential of layer 1 in the lower part of ZH1I and layer 2 in the upper part of ZH1II is moderate. The second layer in the lower part of ZH1II has the poorest properties and lower exploration potential. The main factors affecting the favorable area include physical properties, oil and gas display, sand body thickness, etc.

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.