Large-scale Stochastic Simulations Research Articles

Assessment of long-term trends in genetic mean and variance after the introduction of genomic selection in layers: a simulation study.

Nucleus-based breeding programs are characterized by intense selection that results in high genetic gain, which inevitably means reduction of genetic variation in the breeding population. Therefore, genetic variation in such breeding systems is typically managed systematically, for example, by avoiding mating the closest relatives to limit progeny inbreeding. However, intense selection requires maximum effort to make such breeding programs sustainable in the long-term. The objective of this study was to use simulation to evaluate the long-term impact of genomic selection on genetic mean and variance in an intense layer chicken breeding program. We developed a large-scale stochastic simulation of an intense layer chicken breeding program to compare conventional truncation selection to genomic truncation selection optimized with either minimization of progeny inbreeding or full-scale optimal contribution selection. We compared the programs in terms of genetic mean, genic variance, conversion efficiency, rate of inbreeding, effective population size, and accuracy of selection. Our results confirmed that genomic truncation selection has immediate benefits compared to conventional truncation selection in all specified metrics. A simple minimization of progeny inbreeding after genomic truncation selection did not provide any significant improvements. Optimal contribution selection was successful in having better conversion efficiency and effective population size compared to genomic truncation selection, but it must be fine-tuned for balance between loss of genetic variance and genetic gain. In our simulation, we measured this balance using trigonometric penalty degrees between truncation selection and a balanced solution and concluded that the best results were between 45° and 65°. This balance is specific to the breeding program and depends on how much immediate genetic gain a breeding program may risk vs. save for the future. Furthermore, our results show that the persistence of accuracy is better with optimal contribution selection compared to truncation selection. In general, our results show that optimal contribution selection can ensure long-term success in intensive breeding programs using genomic selection.

Infectious Disease Control in Metapopulations with Limited Resources

Motivated by unique challenges faced in containing the 2014 Ebola outbreak in West Africa, we develop a framework to dynamically allocate limited resources to several possibly connected populations where the disease transmission is stochastic. We formulate this problem as a stochastic dynamic program. However, as the state and action spaces grow exponentially with the size of the problem, the standard solution techniques do not apply. We propose two solution methodologies along with several benchmark policies. The first approach considers a dynamic one-step look-ahead policy which is equivalent to a nonlinear integer knapsack that scales well with the problem size. The second approach is a modification of a myopic incidence policy found in the literature. In addition to testing the proposed policies in a simulation setting of the optimization framework, we develop a large-scale stochastic simulation for 2014 Ebola outbreak in a case study. We calibrate and validate the stochastic simulation model with real-world data from Sierra Leone. Our results provide insights on efficient prioritization and resource allocation in this setting.

Diversity of synaptic protein complexes as a function of the abundance of their constituent proteins: A modeling approach.

The postsynaptic density (PSD) is a dense protein network playing a key role in information processing during learning and memory, and is also indicated in a number of neurological disorders. Efforts to characterize its detailed molecular organization are encumbered by the large variability of the abundance of its constituent proteins both spatially, in different brain areas, and temporally, during development, circadian rhythm, and also in response to various stimuli. In this study we ran large-scale stochastic simulations of protein binding events to predict the presence and distribution of PSD complexes. We simulated the interactions of seven major PSD proteins (NMDAR, AMPAR, PSD-95, SynGAP, GKAP, Shank3, Homer1) based on previously published, experimentally determined protein abundance data from 22 different brain areas and 42 patients (altogether 524 different simulations). Our results demonstrate that the relative ratio of the emerging protein complexes can be sensitive to even subtle changes in protein abundances and thus explicit simulations are invaluable to understand the relationships between protein availability and complex formation. Our observations are compatible with a scenario where larger supercomplexes are formed from available smaller binary and ternary associations of PSD proteins. Specifically, Homer1 and Shank3 self-association reactions substantially promote the emergence of very large protein complexes. The described simulations represent a first approximation to assess PSD complex abundance, and as such, use significant simplifications. Therefore, their direct biological relevance might be limited but we believe that the major qualitative findings can contribute to the understanding of the molecular features of the postsynapse.

Robustness of the microtubule network self-organization in epithelia.

Robustness of biological systems is crucial for their survival, however, for many systems its origin is an open question. Here, we analyze one subcellular level system, the microtubule cytoskeleton. Microtubules self-organize into a network, along which cellular components are delivered to their biologically relevant locations. While the dynamics of individual microtubules is sensitive to the organism's environment and genetics, a similar sensitivity of the overall network would result in pathologies. Our large-scale stochastic simulations show that the self-organization of microtubule networks is robust in a wide parameter range in individual cells. We confirm this robustness in vivo on the tissue-scale using genetic manipulations of Drosophila epithelial cells. Finally, our minimal mathematical model shows that the origin of robustness is the separation of time-scales in microtubule dynamics rates. Altogether, we demonstrate that the tissue-scale self-organization of a microtubule network depends only on cell geometry and the distribution of the microtubule minus-ends.

Meta-model based stochastic simulation of fire barrier cold-side temperature

Meta-modeling helps to reduce the computational cost of large-scale stochastic simulations. It, however, brings froth additional uncertainty in the prediction. In this work, we illustrate the use of the Response Surface Method (RSM) and Gaussian Process (GP) regression to approximate the cold-side temperature of a fire barrier in a stochastic compartment fire scenario. The deterministic model is a 1D Finite Element Method (FEM) solving heat-diffusion. We also illustrate a method for the quantification and compensation of the meta-modeling uncertainty. Additionally, we present a simple technique to create a RSM model with an arbitrary order polynomial function. The results show that both RSM and GP model can produce accurate stochastic simulations despite only few deterministic data points. The simple polynomial-based RSM approximation, however, fails when the heat-transfer is affected by exothermic reactions. This is in contrast with GP, where the combination of various kernel functions makes the prediction flexible enough even for the highly non-linear problem. The meta-modeling uncertainty can be quantified during the testing of meta-models and later be used to improve the accuracy of stochastic predictions.

Assessing leakage detectability at geologic CO2 sequestration sites using the probabilistic collocation method

We present an efficient methodology for assessing leakage detectability at geologic carbon sequestration sites under parameter uncertainty. Uncertainty quantification (UQ) and risk assessment are integral and, in many countries, mandatory components of geologic carbon sequestration projects. A primary goal of risk assessment is to evaluate leakage potential from anthropogenic and natural features, which constitute one of the greatest threats to the integrity of carbon sequestration repositories. The backbone of our detectability assessment framework is the probability collocation method (PCM), an efficient, nonintrusive, uncertainty-quantification technique that can enable large-scale stochastic simulations that are based on results from only a small number of forward-model runs. The metric for detectability is expressed through an extended signal-to-noise ratio (SNR), which incorporates epistemic uncertainty associated with both reservoir and aquifer parameters. The spatially heterogeneous aquifer hydraulic conductivity is parameterized using Karhunen–Loève (KL) expansion. Our methodology is demonstrated numerically for generating probability maps of pressure anomalies and for calculating SNRs. Results indicate that the likelihood of detecting anomalies depends on the level of uncertainty and location of monitoring wells. A monitoring well located close to leaky locations may not always yield the strongest signal of leakage when the level of uncertainty is high. Therefore, our results highlight the need for closed-loop site characterization, monitoring network design, and leakage source detection.

Mutual repression enhances the steepness and precision of gene expression boundaries.

Embryonic development is driven by spatial patterns of gene expression that determine the fate of each cell in the embryo. While gene expression is often highly erratic, embryonic development is usually exceedingly precise. In particular, gene expression boundaries are robust not only against intra-embryonic fluctuations such as noise in gene expression and protein diffusion, but also against embryo-to-embryo variations in the morphogen gradients, which provide positional information to the differentiating cells. How development is robust against intra- and inter-embryonic variations is not understood. A common motif in the gene regulation networks that control embryonic development is mutual repression between pairs of genes. To assess the role of mutual repression in the robust formation of gene expression patterns, we have performed large-scale stochastic simulations of a minimal model of two mutually repressing gap genes in Drosophila, hunchback (hb) and knirps (kni). Our model includes not only mutual repression between hb and kni, but also the stochastic and cooperative activation of hb by the anterior morphogen Bicoid (Bcd) and of kni by the posterior morphogen Caudal (Cad), as well as the diffusion of Hb and Kni between neighboring nuclei. Our analysis reveals that mutual repression can markedly increase the steepness and precision of the gap gene expression boundaries. In contrast to other mechanisms such as spatial averaging and cooperative gene activation, mutual repression thus allows for gene-expression boundaries that are both steep and precise. Moreover, mutual repression dramatically enhances their robustness against embryo-to-embryo variations in the morphogen levels. Finally, our simulations reveal that diffusion of the gap proteins plays a critical role not only in reducing the width of the gap gene expression boundaries via the mechanism of spatial averaging, but also in repairing patterning errors that could arise because of the bistability induced by mutual repression.

A method for the prediction of multitrait breeding values for use in stochastic simulation to compare progeny‐testing schemes, with large progeny groups for proven sires

A method of approximating estimated breeding values (EBV) from a multivariate distribution of true breeding values (TBV) and EBV is proposed for use in large-scale stochastic simulation of alternative breeding schemes with a complex breeding goal. The covariance matrix of the multivariate distributions includes the additive genetic (co)variances and approximated prediction error (co)variances at different selection stages in the life of the animal. The prediction error (co)variance matrix is set up for one animal at a time, utilizing information on the selection candidate and its offspring, the parents, as well as paternal and maternal half- sibs. The EBV are a regression on TBV taking individual uncertainty into account, but with additional 'free' variation drawn at random. With the current information included in the calculation of the prediction error variance of a selection candidate, it is concluded that the method can be used to optimize progeny-testing schemes, where the progeny-tested sires are utilized with large progeny groups, e.g. through artificial insemination.

Predicting the Injury Severity of Pedestrians in Frontal Vehicle Crashes based on Empirical, In-depth Accident Data

For proposed pedestrian protection systems, evaluation of safety benefits is required as an integral part of the design and optimization phases. Stochastic (“Monte-Carlo”) simulation techniques are currently being utilized to predict safety benefits in terms of physics; however, converting physics to human benefits requires injury and fatality risk models. To this end, multivariate predictive models for pedestrian fatalities and for injury severity on the ISS scale are estimated using the US Pedestrian Crash Data Study (PCDS). In addition to collision speed, which is the most important single explanatory variable, age, pedestrian physiological characteristics and vehicle parameters are significant multivariate predictors. The in-sample as well as out-of-sample predictive quality is remarkably high. The models are intended to provide an interface to large-scale stochastic simulation and virtual testing of proposed vehicle-based active safety systems for pedestrian protection.

Accounting for Parameter Uncertainty in Large-Scale Stochastic Simulations with Correlated Inputs

This paper considers large-scale stochastic simulations with correlated inputs having normal-to-anything (NORTA) distributions with arbitrary continuous marginal distributions. Examples of correlated inputs include processing times of workpieces across several workcenters in manufacturing facilities and product demands and exchange rates in global supply chains. Our goal is to obtain mean performance measures and confidence intervals for simulations with such correlated inputs by accounting for the uncertainty around the NORTA distribution parameters estimated from finite historical input data. This type of uncertainty is known as the parameter uncertainty in the discrete-event stochastic simulation literature. We demonstrate how to capture parameter uncertainty with a Bayesian model that uses Sklar's marginal-copula representation and Cooke's copula-vine specification for sampling the parameters of the NORTA distribution. The development of such a Bayesian model well suited for handling many correlated inputs is the primary contribution of this paper. We incorporate the Bayesian model into the simulation replication algorithm for the joint representation of stochastic uncertainty and parameter uncertainty in the mean performance estimate and the confidence interval. We show that our model improves both the consistency of the mean line-item fill-rate estimates and the coverage of the confidence intervals in multiproduct inventory simulations with correlated demands.

Contributors

Contributors

In This Issue

In This Issue

Multiscale modeling through MUSIC multi-simulation: Modeling a dendritic spine using MOOSE and NeuroRD

The nervous system encompasses structure and phenomena at different spatial and temporal scales from molecule to behavior. In addition, different scales are described by different physical and mathematical formalisms. The dynamics of second messenger pathways can be formulated as stochastic reaction-diffusion systems [1] while the electrical dynamics of the neuronal membrane is often described by compartment models and the Hodgkin-Huxley formalism. In neuroscience, there is an increasing need and interest to study multi-scale phenomena where multiple scales and physical formalisms are covered by a single model. While there exists simulators/frameworks, such as GENESIS and MOOSE [3], which span such scales (kinetikit/HH-models), most software applications are specialized for a given domain. Here, we report about initial steps towards a framework for multi-scale modeling which builds on the concept of multi-simulations [2]. We aim to provide a standard API and communication framework allowing parallel simulators targeted at different scales and/or different physics to communicate on-line in a cluster environment. Specifically, we show prototype work on simulating the effect on receptor induced cascades on membrane excitability. Electrical properties of a compartment model is simulated in MOOSE, while receptor induced cascades are simulated in NeuroRD [4,7] . In a prototype system, the two simulators are connected using PyMOOSE [5] and JPype [6]. The two models with their different physical properties (membrane currents in MOOSE, molecular biophysics in NeuroRD), are joined into a single model system. We demonstrate the interaction of the numerical solvers of two simulators (MOOSE, NeuroRD) targeted at different spatiotemporal scales and different physics while solving a multi-scale problem. We analyze some of the problems that may arise in multi-scale multi-simulations and present requirements for a generic API for parallel solvers. This work represents initial steps towards a flexible modular framework for simulation of large-scale multi-scale multi-physics problems in neuroscience.

The Role of Type 4 Phosphodiesterases in Generating Microdomains of cAMP: Large Scale Stochastic Simulations

Cyclic AMP (cAMP) and its main effector Protein Kinase A (PKA) are critical for several aspects of neuronal function including synaptic plasticity. Specificity of synaptic plasticity requires that cAMP activates PKA in a highly localized manner despite the speed with which cAMP diffuses. Two mechanisms have been proposed to produce localized elevations in cAMP, known as microdomains: impeded diffusion, and high phosphodiesterase (PDE) activity. This paper investigates the mechanism of localized cAMP signaling using a computational model of the biochemical network in the HEK293 cell, which is a subset of pathways involved in PKA-dependent synaptic plasticity. This biochemical network includes cAMP production, PKA activation, and cAMP degradation by PDE activity. The model is implemented in NeuroRD: novel, computationally efficient, stochastic reaction-diffusion software, and is constrained by intracellular cAMP dynamics that were determined experimentally by real-time imaging using an Epac-based FRET sensor (H30). The model reproduces the high concentration cAMP microdomain in the submembrane region, distinct from the lower concentration of cAMP in the cytosol. Simulations further demonstrate that generation of the cAMP microdomain requires a pool of PDE4D anchored in the cytosol and also requires PKA-mediated phosphorylation of PDE4D which increases its activity. The microdomain does not require impeded diffusion of cAMP, confirming that barriers are not required for microdomains. The simulations reported here further demonstrate the utility of the new stochastic reaction-diffusion algorithm for exploring signaling pathways in spatially complex structures such as neurons.

Copula-Based Multivariate Input Models for Stochastic Simulation

As large-scale discrete-event stochastic simulation becomes a tool that is used routinely for the design and analysis of stochastic systems, the need for input-modeling support with the ability to represent complex interactions and interdependencies among the components of multivariate time-series input processes is more critical than ever. Motivated by the failure of independent and identically distributed random variables to represent such input processes, a comprehensive framework called Vector-Autoregressive-To-Anything (VARTA) has been introduced for multivariate time-series input modeling. Despite its flexibility in capturing a wide variety of distributional shapes, we show that VARTA falls short in representing dependence structures that arise in situations where extreme component realizations occur together. We demonstrate that it is possible to extend VARTA to work for such dependence structures via the use of the copula theory, which has been used primarily for random vectors in the simulation input-modeling literature, for multivariate time-series input modeling. We show that our copula-based multivariate time-series input model, which includes VARTA as a special case, allows the development of statistically valid fitting and fast sampling algorithms well suited for driving large-scale stochastic simulations.

Controlling endemic cholera with oral vaccines.

BackgroundAlthough advances in rehydration therapy have made cholera a treatable disease with low case-fatality in settings with appropriate medical care, cholera continues to impose considerable mortality in the world's most impoverished populations. Internationally licensed, killed whole-cell based oral cholera vaccines (OCVs) have been available for over a decade, but have not been used for the control of cholera. Recently, these vaccines were shown to confer significant levels of herd protection, suggesting that the protective potential of these vaccines has been underestimated and that these vaccines may be highly effective in cholera control when deployed in mass immunization programs. We used a large-scale stochastic simulation model to investigate the possibility of controlling endemic cholera with OCVs.Methods and FindingsWe construct a large-scale, stochastic cholera transmission model of Matlab, Bangladesh. We find that cholera transmission could be controlled in endemic areas with 50% coverage with OCVs. At this level of coverage, the model predicts that there would be an 89% (95% confidence interval [CI] 72%–98%) reduction in cholera cases among the unvaccinated, and a 93% (95% CI 82%–99%) reduction overall in the entire population. Even a more modest coverage of 30% would result in a 76% (95% CI 44%–95%) reduction in cholera incidence for the population area covered. For populations that have less natural immunity than the population of Matlab, 70% coverage would probably be necessary for cholera control, i.e., an annual incidence rate of ≤ 1 case per 1,000 people in the population.ConclusionsEndemic cholera could be reduced to an annual incidence rate of ≤ 1 case per 1,000 people in endemic areas with biennial vaccination with OCVs if coverage could reach 50%–70% depending on the level of prior immunity in the population. These vaccination efforts could be targeted with careful use of ecological data.

Mitigation strategies for pandemic influenza in the United States

Recent human deaths due to infection by highly pathogenic (H5N1) avian influenza A virus have raised the specter of a devastating pandemic like that of 1917-1918, should this avian virus evolve to become readily transmissible among humans. We introduce and use a large-scale stochastic simulation model to investigate the spread of a pandemic strain of influenza virus through the U.S. population of 281 million individuals for R(0) (the basic reproductive number) from 1.6 to 2.4. We model the impact that a variety of levels and combinations of influenza antiviral agents, vaccines, and modified social mobility (including school closure and travel restrictions) have on the timing and magnitude of this spread. Our simulations demonstrate that, in a highly mobile population, restricting travel after an outbreak is detected is likely to delay slightly the time course of the outbreak without impacting the eventual number ill. For R(0) < 1.9, our model suggests that the rapid production and distribution of vaccines, even if poorly matched to circulating strains, could significantly slow disease spread and limit the number ill to <10% of the population, particularly if children are preferentially vaccinated. Alternatively, the aggressive deployment of several million courses of influenza antiviral agents in a targeted prophylaxis strategy may contain a nascent outbreak with low R(0), provided adequate contact tracing and distribution capacities exist. For higher R(0), we predict that multiple strategies in combination (involving both social and medical interventions) will be required to achieve similar limits on illness rates.

A novel approach to model neuronal signal transduction using stochastic differential equations

We introduce a new approach to model the behavior of neuronal signal transduction networks using stochastic differential equations. We present first a mathematical formulation for a stochastic model of protein kinase C pathway. Different kinds of numerical integration methods, including the explicit and implicit Euler–Maruyama methods, are used to solve the It o ^ form of the stochastic model. Stochastic models may provide more realistic representations for the chemical species in signal transduction networks compared to deterministic models. Our approach has the advantage of being computationally less demanding in the context of large-scale stochastic simulations than other approaches where individual chemical interactions are simulated stochastically.

Epitaxial mounding in limited-mobility models of surface growth

We study, through large-scale stochastic simulations using the noise reduction technique, surface growth via vapor deposition, e.g., molecular beam epitaxy (MBE), for simple nonequilibrium limited mobility solid-on-solid growth models, such as the Family model, the Das Sarma--Tamborenea model, the Wolf-Villain (WV) model, the larger-curvature (LC) model, and other related models. We find that $(d=2+1)$-dimensional surface growth in several noise reduced models (most notably the WV and LC models) exhibits spectacular quasiregular mound formation with slope selection in their dynamical surface morphology in contrast to the standard statistically scale-invariant kinetically rough surface growth expected (and earlier reported in the literature) for such growth models. The mounding instability in these epitaxial growth models does not involve the Ehrlich-Schwoebel step-edge diffusion barrier. The mounded morphology in these growth models arises from the interplay between the line tension along step edges in the plane parallel to the average surface and the suppression of noise and island nucleation. The line tension tends to stabilize some of the step orientations that coincide with in-plane high-symmetry crystalline directions, and thus the mounds that are formed assume quasiregular structures. The noise reduction technique developed originally for Eden-type models can be used to control the stochastic noise and enhance diffusion along the step edge, which ultimately leads to the formation of quasiregular mounds during growth. We show that by increasing the diffusion surface length together with supression of nucleation and deposition noise, one can obtain a self-organization of the pyramids in quasiregular patterns. The mounding instability in these simple epitaxial growth models is closely related to the cluster-edge diffusion (as opposed to step-edge barrier) driven mounding in MBE growth, which has been recently discussed in the literature. The epitaxial mound formation studied here is a kinetic-topological instability [which can happen only in $(d=2+1)$-dimensional, or higher dimensional, growth, but not in $(d=1+1)$-dimensional growth because no cluster diffusion around a closed surface loop is possible in ``one-dimensional'' surfaces], which is likely to be quite generic in real MBE-type surface growth. Our extensive numerical simulations produce mounded (and slope-selected) surface growth morphologies which are strikingly visually similar to many recently reported experimental MBE growth morphologies.