Stochastic Spatial Simulations Research Articles

Structure of the HIV immature lattice allows for essential lattice remodeling within budded virions.

For HIV virions to become infectious, the immature lattice of Gag polyproteins attached to the virion membrane must be cleaved. Cleavage cannot initiate without the protease formed by the homo-dimerization of domains linked to Gag. However, only 5% of the Gag polyproteins, termed Gag-Pol, carry this protease domain, and they are embedded within the structured lattice. The mechanism of Gag-Pol dimerization is unknown. Here, we use spatial stochastic computer simulations of the immature Gag lattice as derived from experimental structures, showing that dynamics of the lattice on the membrane is unavoidable due to the missing 1/3 of the spherical protein coat. These dynamics allow for Gag-Pol molecules carrying the protease domains to detach and reattach at new places within the lattice. Surprisingly, dimerization timescales of minutes or less are achievable for realistic binding energies and rates despite retaining most of the large-scale lattice structure. We derive a formula allowing extrapolation of timescales as a function of interaction free energy and binding rate, thus predicting how additional stabilization of the lattice would impact dimerization times. We further show that during assembly, dimerization of Gag-Pol is highly likely and therefore must be actively suppressed to prevent early activation. By direct comparison to recent biochemical measurements within budded virions, we find that only moderately stable hexamer contacts (-12kBT<∆G<-8kBT) retain both the dynamics and lattice structures that are consistent with experiment. These dynamics are likely essential for proper maturation, and our models quantify and predict lattice dynamics and protease dimerization timescales that define a key step in understanding formation of infectious viruses.

Role of predator-prey reversal in rock-paper-scissors models

In this letter we consider a single parameter generalization of the standard three species Rock-Paper-Scissors (RPS) model allowing for predator-prey reversal. This model, which shall be referred to as κRPS model, incorporates bidirectional predator-prey interactions between all the species in addition to the unidirectional predator-prey interactions of the standard RPS model. We study the dynamics of a May-Leonard formulation of the κRPS model using lattice-based spatial stochastic simulations with random initial conditions. We find that if the simulation lattices are sufficiently large for the coexistence of all three species to be maintained, the model asymptotically leads to the formation of spiral patterns whose evolution is qualitatively similar to that of the standard RPS model, albeit with larger characteristic length and time scales. We show that if the likelihood of predator-prey reversal is sufficiently large there are two distinct scaling regimes: one transient curvature dominated regime in which the characteristic length of the population network grows with time and another where it becomes a constant. We also estimate the dependence of the asymptotic value of the characteristic length of the population network on the likelihood of predator-prey reversal and show that if the simulation lattices are not sufficiently large then predator-prey reversal could potentially have a negative impact on coexistence. Finally, we interpret these results by considering the much simpler dynamics of circular domains.

Stochastic Spatial Binary Simulation with Multivariate Normal Distribution for Illustrating Future Evolution of Umbrella-Shape Summer Shelter under Climate Change

Surface temperature has increased due to the impact of climate change, and the related weather events, such as heat waves and urban heat island, are occurring more frequently than before. Local governments and planners consider these impacts of climate change and try to avoid them. One of the mainly used structural tools is building summer shelters. A critical issue for urban planners is to test how many shelters should be added and how to distribute the shelters to cope with the impact of climate change. Stochastic simulation models can be a good option to randomize locations of shelters and to see how beneficial for living the shelters can be. Therefore, a novel stochastic simulation model is proposed for distributing summer shelters for coping with the climate change impact. This study proposes a stochastic spatial binary simulation with multivariate normal distribution (SSBM) which contains two major procedures consisting of (1) simulation-based derivation of the empirical function and (2) stochastic simulation of spatial binary data with multivariate normal distribution and the derived empirical function. The proposed model is applied to a case study in Jinju City, South Korea, for the umbrella-shape summer shelters (USS). Results concluded that the proposed SSBM reproduced the statistical characteristics of USS and can be a good alternative to model the locations of USS, including the impact of climate change and investigating the evolution of the USS in the future.

How local antipredator response unbalances the rock-paper-scissors model.

Antipredator behavior is a self-preservation strategy present in many biological systems, where individuals join the effort in a collective reaction to avoid being caught by an approaching predator. We study a nonhierarchical tritrophic system, whose predator-prey interactions are described by the rock-paper-scissors game rules. We perform a set of spatial stochastic simulations where organisms of one out of the species can resist predation in a collective strategy. The drop in predation capacity is local, which means that each predator faces a particular opposition depending on the prey group size surrounding it. Considering that the interference in a predator action depends on the prey's physical and cognitive ability, we explore the role of a conditioning factor that indicates the fraction of the species apt to perform the antipredator strategy. Because of the local unbalancing of the cyclic predator-prey interactions, departed spatial domains mainly occupied by a single species emerge. Unlike the rock-paper-scissors model with a weak species because of a nonlocal reason, our findings show that if the predation probability of one species is reduced because individuals face local antipredator response, the species does not predominate. Instead, the local unbalancing of the rock-paper-scissors model results in the prevalence of the weak species' prey. Finally, the outcomes show that local unevenness may jeopardize biodiversity, with the coexistence being more threatened for high mobility.

Did Maxwell Dream of Electrical Bacteria?

We propose a model for bacterial Quorum Sensing based on an auxiliary electrostatic-like interaction originating from a fictitious electrical charge that represents bacteria activity. A cooperative mechanism for charge/activity exchange is introduced to implement chemotaxis and replication. The bacteria system is thus represented by means of a complex resistor network where link resistances take into account the allowed activity-flow among individuals. By explicit spatial stochastic simulations, we show that the model exhibits different quasi-realistic behaviors from colony formation to biofilm aggregation. The electrical signal associated with Quorum Sensing is analyzed in space and time and provides useful information about the colony dynamics. In particular, we analyze the transition between the planktonic and colony phases as the intensity of Quorum Sensing is varied.

Parity effects in rock-paper-scissors type models with a number of species [formula omitted

We investigate the impact of parity on the abundance of weak species in the context of the simplest generalization of the rock-paper-scissors model to an arbitrary number of species — we consider models with a total number of species (NS) between 3 and 12, having one or more (weak) species characterized by a reduced predation probability (by a factor of Pw with respect to the other species). We show, using lattice based spatial stochastic simulations with random initial conditions, large enough for coexistence to prevail, that parity effects are significant. We find that the performance of weak species is dependent on whether the total number of species is even or odd, especially for NS≤8, with odd numbers of species being on average more favourable to weak species than even ones. We further show that, despite the significant dispersion observed among individual models, a weak species has on average a higher abundance than a strong one if Pw is sufficiently smaller than unity — the notable exception being the four species case.

Mobility-limiting antipredator response in the rock-paper-scissors model

Antipredator behavior is present in many biological systems where individuals collectively react to an imminent attack. The antipredator response may influence spatial pattern formation and ecosystem stability but requires an organism's cost to contribute to the collective effort. We investigate a nonhierarchical tritrophic system, whose predator-prey interactions are described by the rock-paper-scissors game rules. In our spatial stochastic simulations, the radius of antipredator response defines the maximum prey group size that disturbs the predator's action, determining the individual cost to participate in antipredator strategies. We consider that each organism contributes equally to the collective effort, having its mobility limited by the proportion of energy devoted to the antipredator reaction. Our outcomes show that the antipredator response leads to spiral patterns, with the segregation of organisms of the same species occupying departed spatial domains. We found that a less localized antipredator response increases the average size of the single-species patches, improving the protection of individuals against predation. Finally, our findings show that although the increase of the predation risk for a more localized antipredator response, the high mobility constraining benefits species coexistence. Our results may help ecologists understand the mechanisms leading to the stability of biological systems where locality is crucial to behavioral interactions among species.

What Are the Human Resources Required to Control a Foot-and-Mouth Disease Outbreak in Austria?

Contingency planning allows veterinary authorities to prepare a rapid response in the event of a disease outbreak. A recently published foot-and-mouth disease (FMD) simulation study indicated concerns whether capacity was sufficient to control a potential FMD epidemic in Austria. The objectives of the study presented here were to estimate the human resources required to implement FMD control measures and to identify areas of the operational activities that could potentially delay successful control of the disease. The stochastic spatial simulation model EuFMDiS (The European Foot-and-Mouth Disease Spread Model) was used to simulate a potential FMD outbreak and its economic impact, including different control scenarios based on variations of culling, vaccination, and pre-emptive depopulation. In this context, the utilization of human resources was assessed based on the associated EuFMDiS output regarding the performance of operational activities. The assessments show that the number of personnel needed in an outbreak with a stamping-out policy would reach the peak at the end of the second week of control with a median of 540 (257–926) individuals, out of which 31% would be veterinarians. Approximately 58% of these human resources would be attributable to surveillance, followed by staff for cleaning and disinfection activities. Our analysis demonstrates that, of the operational activities, surveillance personnel were the largest factor influencing the magnitude of the outbreak. The aim of the assessment presented here is to assist veterinary authorities in the contingency planning of required human resources to respond effectively to an outbreak of animal diseases such as FMD.

Spatially variable pesticide application in olive groves: Evaluation of potential pesticide-savings through stochastic spatial simulation algorithms.

Site-specific management using spatial crown volume characterization can greatly reduce the amount of pesticides applied in agricultural treatments performed with air-assisted sprayers, while helping farmers achieve the European legislation on safe use of pesticides. Nevertheless, variable rate treatments in olive groves have received little attention. Thus, field research was conducted in a 20.6-ha traditional olive grove. Two attributes of the trees - tree crown volume (V) and tree projected area - were determined, using 67 samples for V and all trees of the field (1433) for tree projected area. Spatial continuity of both attributes was modelled with exponential variograms. To gain a measure of local uncertainty, stochastic simulation algorithms were applied. One hundred simulated images were obtained for tree projected area using direct sequential simulation. Tree projected area simulations were used to improve spatial prediction of V, more difficult and more expensive to obtain, taking advantage of the high linear correlation between both variables (rxy=0.72,p<0.001). Thus, direct sequential cosimulation was employed to predict the spatial distribution of V, obtaining 100 geostatistical realizations of V. In order to estimate the potential reduction of pesticide use in the farm with variable rate treatments, two cut-off values of V were considered (50 and 100m3crown volume). Local uncertainty, understood as the probability of each tree belonging to a given crown volume interval was determined. Probability maps were further transformed to morphological maps and finally to variable prescription maps. Two scenarios with 2 and 3 management zones (MZs) were obtained. In comparison with a conventional phytosanitary application, the variable rate treatments could reduce the pesticide amounts by 21.3% with 2 MZs, and by 38% with 3 MZs. The joint use of V and tree projected area in stochastic sequential simulation algorithms has shown to be useful to determine MZs in olive groves.

Stochastic Spatial Simulation of Genetic Information Processes in the Minimal Cell

JCVI-syn3A is a genetically minimal cell whose genome consists of 493 genes, of which 452 are protein-coding. The genetic information processing reactions in JCVI-syn3A were recently solved stochastically using a well-stirred assumption which predicts the times of DNA replication initiation and elongation to be 8 and 50 minutes on average respectively after averaging over 1,000 replicate cells. This model approximately doubles the number of proteins according to genome-wide proteomics data as well as the number of ribosomal components over the course of a cell cycle. We introduce spatial heterogeneity and diffusion by constructing a 3D spatial model of JCVI-syn3A using cryo-electron tomograms to infer placement of large cellular components such as membrane boundaries and ribosomes. The ensemble of DNA configurations for the circular chromosome were generated using a random polymer where self-intersections and intersections with ribosomes were forbidden. Membrane-bound degradosome complexes were randomly distributed across the interior of the cell membrane. We fixed positions of the DNA, degradosomes, and ribosomes while simulating diffusion of RNA and proteins with more than 3,000 reactions including transcription of all 452 protein-coding genes and translation and degradation of all 452 mRNA under a reaction-diffusion master equation formulation using our GPU-based software, Lattice Microbes. Introducing spatial heterogeneity and requiring mRNA to diffuse to ribosomes and degradosomes to be translated or decayed results in a distribution of mRNA half-lives in agreement with measurements in other gram-positive bacteria Mycoplasma gallisepticum and Bacillus subtilis.

Performance of weak species in the simplest generalization of the rock-paper-scissors model to four species

We investigate the problem of the predominance and survival of "weak" species in the context of the simplest generalization of the spatial stochastic rock-paper-scissors model to four species by considering models in which one, two, or three species have a reduced predation probability. We show, using lattice based spatial stochastic simulations with random initial conditions, that if only one of the four species has its probability reduced, then the most abundant species is the prey of the "weakest" (assuming that the simulations are large enough for coexistence to prevail). Also, among the remaining cases, we present examples in which "weak" and "strong" species have similar average abundances and others in which either of them dominates-the most abundant species being always a prey of a weak species with which it maintains a unidirectional predator-prey interaction. However, in contrast to the three-species model, we find no systematic difference in the global performance of weak and strong species, and we conjecture that a similar result will hold if the number of species is further increased. We also determine the probability of single species survival and coexistence as a function of the lattice size, discussing its dependence on initial conditions and on the change to the dynamics of the model which results from the extinction of one of the species.

Translating surveillance data into incidence estimates.

Monitoring a population for a disease requires the hosts to be sampled and tested for the pathogen. This results in sampling series from which we may estimate the disease incidence, i.e. the proportion of hosts infected. Existing estimation methods assume that disease incidence does not change between monitoring rounds, resulting in an underestimation of the disease incidence. In this paper, we develop an incidence estimation model accounting for epidemic growth with monitoring rounds that sample varying incidence. We also show how to accommodate the asymptomatic period that is the characteristic of most diseases. For practical use, we produce an approximation of the model, which is subsequently shown to be accurate for relevant epidemic and sampling parameters. Both the approximation and the full model are applied to stochastic spatial simulations of epidemics. The results prove their consistency for a very wide range of situations. The estimation model is made available as an online application. This article is part of the theme issue 'Modelling infectious disease outbreaks in humans, animals and plants: epidemic forecasting and control'. This theme issue is linked with the earlier issue 'Modelling infectious disease outbreaks in humans, animals and plants: approaches and important themes'.

EGFRD in all dimensions.

Biochemical reactions often occur at low copy numbers but at once in crowded and diverse environments. Space and stochasticity therefore play an essential role in biochemical networks. Spatial-stochastic simulations have become a prominent tool for understanding how stochasticity at the microscopic level influences the macroscopic behavior of such systems. While particle-based models guarantee the level of detail necessary to accurately describe the microscopic dynamics at very low copy numbers, the algorithms used to simulate them typically imply trade-offs between computational efficiency and biochemical accuracy. eGFRD (enhanced Green's Function Reaction Dynamics) is an exact algorithm that evades such trade-offs by partitioning the N-particle system into M ≤ N analytically tractable one- and two-particle systems; the analytical solutions (Green's functions) then are used to implement an event-driven particle-based scheme that allows particles to make large jumps in time and space while retaining access to their state variables at arbitrary simulation times. Here we present "eGFRD2," a new eGFRD version that implements the principle of eGFRD in all dimensions, thus enabling efficient particle-based simulation of biochemical reaction-diffusion processes in the 3D cytoplasm, on 2D planes representing membranes, and on 1D elongated cylinders representative of, e.g., cytoskeletal tracks or DNA; in 1D, it also incorporates convective motion used to model active transport. We find that, for low particle densities, eGFRD2 is up to 6 orders of magnitude faster than conventional Brownian dynamics. We exemplify the capabilities of eGFRD2 by simulating an idealized model of Pom1 gradient formation, which involves 3D diffusion, active transport on microtubules, and autophosphorylation on the membrane, confirming recent experimental and theoretical results on this system to hold under genuinely stochastic conditions.

A Stochastic Spatial Simulation Method for Self-Assembly Reactions

The ability of collections of molecules to spontaneously assemble into large functional complexes is central to all cellular processes. Yet the fine details of complex self-assembly processes normally cannot be directly experimentally observed due to their small scale and rapid dynamics. Intractability of experimental approaches is particularly acute for understanding self-assembly in vivo, which may operate quite differently from purified in vitro models due to such effects as spatial confinement, macromolecular crowding, and influences of extrinsic cellular factors. Simulation has thus proven a valuable adjunct to experimental methods, facilitating interpretation of experimental data and inferences about experimentally unobservable aspects of self-assembly. A variety of simulation methodologies have been used for modeling self-assembly dynamics --- including approaches based on Differential Equations (DE), Gillespie Algorithm (SSA), Brownian Dynamics (BD), Green's Function Reaction Dynamics (GFRD), and variations thereon --- each offering tradeoffs with respect to the ranges of phenomena they can model, their computational tractability, and the difficulty of fitting them to experimental measurements. Here, we explore a multiscale approach combining efficient SSA-like sampling suitable for well-mixed systems with the slower but space-aware GFRD model, assuming as with GFRD that reactions occur in a spatially heterogeneous environment that must be explicitly modeled. We generalize the GFRD model to more efficiently handle bimolecular association reactions fundamental to self-assembly, using a rejection method for exact and efficient sampling of reaction locations in spatially heterogeneous regions. To further accelerate models with minimal loss of accuracy, we are exploring regression-based methods for approximate integration. While the resulting model remains a theoretical abstraction, future work will explore its ability to capture an important subset of cellular self-assembly for which classic SSA is too simplified but GFRD and BD too inefficient to capture accurately the phenomena in realistic biological parameter ranges.

Single-cell variability in multicellular organisms

Noisy gene expression is of fundamental importance to single cells, and is therefore widely studied in single-celled organisms. Extending these studies to multicellular organisms is challenging since their cells are generally not isolated, but individuals in a tissue. Cell–cell coupling via signalling, active transport or pure diffusion, ensures that tissue-bound cells are neither fully independent of each other, nor an entirely homogeneous population. In this article, we show that increasing the strength of coupling between cells can either increase or decrease the single-cell variability (and, therefore, the heterogeneity of the tissue), depending on the statistical properties of the underlying genetic network. We confirm these predictions using spatial stochastic simulations of simple genetic networks, and experimental data from animal and plant tissues. The results suggest that cell–cell coupling may be one of several noise-control strategies employed by multicellular organisms, and highlight the need for a deeper understanding of multicellular behaviour.

Compartmental and Spatial Rule-Based Modeling with Virtual Cell

In rule-based modeling, molecular interactions are systematically specified in the form of reaction rules that serve as generators of reactions. This provides a way to account for all the potential molecular complexes and interactions among multivalent or multistate molecules. Recently, we introduced rule-based modeling into the Virtual Cell (VCell) modeling framework, permitting graphical specification of rules and merger of networks generated automatically (using the BioNetGen modeling engine) with hand-specified reaction networks. VCell provides a number of ordinary differential equation and stochastic numerical solvers for single-compartment simulations of the kinetic systems derived from these networks, and agent-based network-free simulation of the rules. In this work, compartmental and spatial modeling of rule-based models has been implemented within VCell. To enable rule-based deterministic and stochastic spatial simulations and network-free agent-based compartmental simulations, the BioNetGen and NFSim engines were each modified to support compartments. In the new rule-based formalism, every reactant and product pattern and every reaction rule are assigned locations. We also introduce the rule-based concept of molecular anchors. This assures that any species that has a molecule anchored to a predefined compartment will remain in this compartment. Importantly, in addition to formulation of compartmental models, this now permits VCell users to seamlessly connect reaction networks derived from rules to explicit geometries to automatically generate a system of reaction-diffusion equations. These may then be simulated using either the VCell partial differential equations deterministic solvers or the Smoldyn stochastic simulator.

Rebinding in biochemical reactions on membranes

The behavior of many biochemical processes depends crucially on molecules rapidly rebinding after dissociating. In the case of multisite protein modification, the importance of rebinding has been demonstrated both experimentally and through several recent computational studies involving stochastic spatial simulations. As rebinding stems from spatio-temporal correlations, theorists have resorted to models that explicitly include space to properly account for the effects of rebinding. However, for reactions in three space dimensions it was recently shown that well-mixed ordinary differential equation (ODE) models can incorporate rebinding by adding connections to the reaction network. The rate constants for these new connections involve the probability that a pair of molecules rapidly rebinds after dissociation. In order to study biochemical reactions on membranes, in this paper we derive an explicit formula for this rebinding probability for reactions in two space dimensions. We show that ODE models can use the formula to replicate detailed stochastic spatial simulations, and that the formula can predict ultrasensitivity for reactions involving multisite modification of membrane-bound proteins. Further, we compute a new concentration-dependent rebinding probability for reactions in three space dimensions. Our analysis predicts that rebinding plays a much larger role in reactions on membranes compared to reactions in cytoplasm.

Parallel STEPS: Large Scale Stochastic Spatial Reaction-Diffusion Simulation with High Performance Computers.

Stochastic, spatial reaction-diffusion simulations have been widely used in systems biology and computational neuroscience. However, the increasing scale and complexity of models and morphologies have exceeded the capacity of any serial implementation. This led to the development of parallel solutions that benefit from the boost in performance of modern supercomputers. In this paper, we describe an MPI-based, parallel operator-splitting implementation for stochastic spatial reaction-diffusion simulations with irregular tetrahedral meshes. The performance of our implementation is first examined and analyzed with simulations of a simple model. We then demonstrate its application to real-world research by simulating the reaction-diffusion components of a published calcium burst model in both Purkinje neuron sub-branch and full dendrite morphologies. Simulation results indicate that our implementation is capable of achieving super-linear speedup for balanced loading simulations with reasonable molecule density and mesh quality. In the best scenario, a parallel simulation with 2,000 processes runs more than 3,600 times faster than its serial SSA counterpart, and achieves more than 20-fold speedup relative to parallel simulation with 100 processes. In a more realistic scenario with dynamic calcium influx and data recording, the parallel simulation with 1,000 processes and no load balancing is still 500 times faster than the conventional serial SSA simulation.

Modeling the Role of Epitope Arrangement on Antibody Binding Stoichiometry in Flaviviruses

Cryo-electron-microscopy (cryo-EM) structures of flaviviruses reveal significant variation in epitope occupancy across different monoclonal antibodies that have largely been attributed to epitope-level differences in conformation or accessibility that affect antibody binding. The consequences of these variations for macroscopic properties such as antibody binding and neutralization are the results of the law of mass action—a stochastic process of innumerable binding and unbinding events between antibodies and the multiple binding sites on the flavivirus in equilibrium—that cannot be directly imputed from structure alone. We carried out coarse-grained spatial stochastic binding simulations for nine flavivirus antibodies with epitopes defined by cryo-EM or x-ray crystallography to assess the role of epitope spatial arrangement on antibody-binding stoichiometry, occupancy, and neutralization. In our simulations, all epitopes were equally competent for binding, representing the upper limit of binding stoichiometry that results from epitope spatial arrangement alone. Surprisingly, our simulations closely reproduced the relative occupancy and binding stoichiometry observed in cryo-EM, without having to account for differences in epitope accessibility or conformation, suggesting that epitope spatial arrangement alone may be sufficient to explain differences in binding occupancy and stoichiometry between antibodies. Furthermore, we found that there was significant heterogeneity in binding configurations even at saturating antibody concentrations, and that bivalent antibody binding may be more common than previously thought. Finally, we propose a structure-based explanation for the stoichiometric threshold model of neutralization.

Simulation of weathered profiles coupled with multivariate block-support simulation of the Puma nickel laterite deposit, Brazil

Modelling and assessing spatial variability and uncertainty of mineral deposits is critical for both capital investments in mining projects as well as operational issues once a mine is developed. However, traditional approaches for modelling geological domains and geostatistical estimation provide smoothed representations of the pertinent deposit attributes, ignore spatial variability and, thereby, can mislead downstream decisions. Spatial variability and related uncertainty in modelling mineral deposit characteristics of interest, ranging from metal content and geological boundaries to geomechanical - geotechnical rock properties, can be modelled and quantified by stochastic spatial simulations. This is demonstrated through a detailed, step- by-step application to the Puma deposit, a major nickel lateritic asset in Brazil, part of the Onça-Puma mining complex. To integrate the variability of the regolith profiles of the deposit, their thicknesses are calculated after an unwrinkling process is applied and the deposit is then jointly simulated using min/max autocorrelation factors (MAF). The realizations serve as geological boundaries within which Ni, Co, Fe, SiO2, MgO and Dry-tonnage factor (DTF) are subsequently jointly simulated directly at block support scale. The final result is a series of equally probable representations of the Puma deposit, which are used to quantify and assess the uncertainty about key aspects of the project at the Puma deposit, such as the uncertainty of the in-situ resources documented herein, and the strict control of the ore's quality that feeds the ferronickel processing plant. The framework presented shows the advantages of the MAF and direct block simulation approaches for the efficient joint simulation of spatially variant geological attributes of large deposits for industrial environments. The methods presented are general, new in the context of geotechnical and geomechanical engineering, and can assist in the modelling of spatial variability and quantification of uncertainty linked to geotechnical rock properties and pertinent lithological boundaries.