Journal articles on the topic 'Gillespie simulations'

Abstract Monte Carlo simulations have been utilized to make a comparative study between diffusion approximation (DA) and the Gillespie algorithm and its dependence on population in the information diffusion model. Diffusion approximation is one of the widely used approximation methods which have been applied in queuing systems, biological systems and other fields. The Gillespie algorithm, on the other hand, is used for simulating stochastic systems. In this article, the validity of diffusion approximation has been studied in relation to the Gillespie algorithm for varying population sizes. It is found that diffusion approximation results in large fluctuations which render forecasting unreliable particularly for a small population. The relative fluctuations in relation to diffusion approximation, as well as to the Gillespie algorithm have been analyzed. To carry out the study, a nonlinear stochastic model of innovation diffusion in a finite population has been considered. The nonlinearity of the problem necessitates use of approximation methods to understand the dynamics of the system. A stochastic differential equation (SDE) has been used to model the innovation diffusion process, and corresponding sample paths have been generated using Monte Carlo simulation methods.

Abstract. The evolution of two-dimensional drop distributions is simulated in this study using a Monte Carlo method.~The stochastic algorithm of Gillespie (1976) for chemical reactions in the formulation proposed by Laurenzi et al. (2002) was used to simulate the kinetic behavior of the drop population. Within this framework species are defined as droplets of specific size and aerosol composition. The performance of the algorithm was checked by comparing the numerical with the analytical solutions found by Lushnikov (1975). Very good agreement was observed between the Monte Carlo simulations and the analytical solution. Simulation results are presented for bi-variate constant and hydrodynamic kernels. The algorithm can be easily extended to incorporate various properties of clouds such as including several crystal habits, different types of soluble CCN, particle charging and drop breakup.

Bacterial heteroresistance (i.e., the co-existence of several subpopulations with different antibiotic susceptibilities) can delay the clearance of bacteria even with long antibiotic exposure. Some proposed mechanisms have been successfully described with mathematical models of drug-target binding where the mechanism’s downstream of drug-target binding are not explicitly modeled and subsumed in an empirical function, connecting target occupancy to antibiotic action. However, with current approaches it is difficult to model mechanisms that involve multi-step reactions that lead to bacterial killing. Here, we have a dual aim: first, to establish pharmacodynamic models that include multi-step reaction pathways, and second, to model heteroresistance and investigate which molecular heterogeneities can lead to delayed bacterial killing. We show that simulations based on Gillespie algorithms, which have been employed to model reaction kinetics for decades, can be useful tools to model antibiotic action via multi-step reactions. We highlight the strengths and weaknesses of current models and Gillespie simulations. Finally, we show that in our models, slight normally distributed variances in the rates of any event leading to bacterial death can (depending on parameter choices) lead to delayed bacterial killing (i.e., heteroresistance). This means that a slowly declining residual bacterial population due to heteroresistance is most likely the default scenario and should be taken into account when planning treatment length.

Plants are vital for man and many species. They are sources of food, medicine, fiber for clothes and materials for shelter. They are a fundamental part of a healthy environment. However, plants are subject to virus diseases. In plants most of the virus propagation is done by a vector. The traditional way of controlling the insects is to use insecticides that have a negative effect on the environment. A more environmentally friendly way to control the insects is to use predators that will prey on the vector, such as birds or bats. In this paper we modify a plant-virus propagation model with delays. The model is written using delay differential equations. However, it can also be expressed in terms of biochemical reactions, which is more realistic for small populations. Since there are always variations in the populations, errors in the measured values and uncertainties, we use two methods to introduce randomness: stochastic differential equations and the Gillespie algorithm. We present numerical simulations. The Gillespie method produces good results for plant-virus population models.

Abstract. The evolution of two-dimensional drop distributions is simulated in this study using a Monte Carlo method. The stochastic algorithm of Gillespie (1976) for chemical reactions in the formulation proposed by Laurenzi et al. (2002) was used to simulate the kinetic behavior of the drop population. Within this framework, species are defined as droplets of specific size and aerosol composition. The performance of the algorithm was checked by a comparison with the analytical solutions found by Lushnikov (1975) and Golovin (1963) and with finite difference solutions of the two-component kinetic collection equation obtained for the Golovin (sum) and hydrodynamic kernels. Very good agreement was observed between the Monte Carlo simulations and the analytical and numerical solutions. A simulation for realistic initial conditions is presented for the hydrodynamic kernel. As expected, the aerosol mass is shifted from small to large particles due to collection process. This algorithm could be extended to incorporate various properties of clouds such several crystals habits, different types of soluble CCN, particle charging and drop breakup.

In this paper, we propose a novel Poisson process generator that uses multiple thermal noise amplifiers (TNAs) as a source of randomness and controls its event rate via a frequency-locked loop (FLL). The increase in the number of TNAs extends the effective bandwidth of amplified thermal noise and hence enhances the maximum event rate the proposed architecture can generate. Verilog-A simulation of the proposed Poisson process generator shows that its maximum event rate can be increased by a factor of 26.5 when the number of TNAs increases from 1 to 10. In order to realize parallel stochastic simulations of the biochemical reaction network, we present a fundamental reaction building block with continuous-time multiplication and addition using an AND gate and a 1-bit current-steering digital-to-analog converter, respectively. Stochastic biochemical reactions consisting of the fundamental reaction building blocks are simulated in Verilog-A, demonstrating that the simulation results are consistent with those of conventional Gillespie algorithm. An increase in the number of TNAs to accelerate the Poisson events and the use of digital AND gates for robust reaction rate calculations allow for faster and more accurate stochastic simulations of biochemical reactions than previous parallel stochastic simulators.

Most recently, the production of polystyrene by Free Radical Polymerization (FRP) via microchannels has been a subject of core interest due to the efficiency of a micro-or milli-reactor brings. In addition, especially in pilot experimentations, a micro or milli-reactor has been known widely to be efficient in monitoring the microstructural end-use features or properties of the polymer as the chain propagates and ultimately terminates. However, the limitations posed by using micro or milli-reactors in process intensification such as clogging of pores can be a bottleneck when tracking the common phenomena associated with FRP such as cage, gel, and glass effects. In this work, the simulation of the synthesis of polystyrene in FRP via microchannels is computed using a robust and time-efficient hybrid Gillespie Algorithm (GA) or Hybrid Stochastic Simulation Algorithm (HSSA). The obtained results of the end-use properties of polystyrene such as Monomer conversion, Polydispersity Index, Number-Average Molar Mass and Weight Average Molar Mass were compared to experimental data. The simulation results agree well with the experimental results reported in this work. Hence, stochastic simulations prove to be an effective tool in making decisions in the context of process intensification of chain growth polymerization reactions even at a large scale.

The elevation of Synthetic Biology from single cells to multicellular simulations would be a significant scale-up. The spatiotemporal behavior of cellular populations has the potential to be prototyped in silico for computer assisted design through ergonomic interfaces. Such a platform would have great practical potential across medicine, industry, research, education and accessible archiving in bioinformatics. Existing Synthetic Biology CAD systems are considered limited regarding population level behavior, and this work explored the in silico challenges posed from biological and computational perspectives. Retaining the connection to Synthetic Biology CAD, an extension of the Infobiotics Workbench Suite was considered, with potential for the integration of genetic regulatory models and/or chemical reaction networks through Next Generation Stochastic Simulator (NGSS) Gillespie algorithms. These were executed using SBML models generated by in-house SBML-Constructor over numerous topologies and benchmarked in association with multicellular simulation layers. Regarding multicellularity, two ground-up multicellular solutions were developed, including the use of Unreal Engine 4 contrasted with CPU multithreading and Blender visualization, resulting in a comparison of real-time versus batch-processed simulations. In conclusion, high-performance computing and client–server architectures could be considered for future works, along with the inclusion of numerous biologically and physically informed features, whilst still pursuing ergonomic solutions.

We discuss how translocation properties of the 20S proteasome influence its length distribution, one of its most important feature for the normal functioning of the immune system. For this we consider a simple one-channel proteasome model and assume that the protein transport depends significantly on the length of a protein located inside the proteasome chamber. Using the master equation approach we show analytically that the length distribution with one dominating peak, observed in the experiments, can be achieved if the transport rate function is in a certain relation with cleavage probabilities and the geometry of a proteasome. Our analytical results are confirmed by numerical simulations of the protein degradation by the proteasome performed using the modified Gillespie algorithm.

Cell spreading plays an important role in the modulation of physiological functions such as inflammation and cancer metastasis. The Brownian ratchet model and Bell's model have been used to simulate actin dynamics and bond kinetics for focal adhesion dynamics, respectively. In the present study, these models were modified and two additional subcellular mechanisms, integrin and myosin kinetics, were incoporated. An integrin recruitment function was introduced to determine the size of a focal adhesion associated with the substrate stiffness. The relationship between myosin concentration and the actin protrusion velocity was described by a first-order differential equation. Subcellular processes, including cell protrusion, focal adhesion formation, and stress fiber formation, were integrated into an axial-symmetric biophysical model, while inputs to the model were kinematic data from time-lapse experiments. Numerical simulations of the model using the Gillespie algorithm showed that dynamics of cell spreading can be well described by the model.

Populations of ecological systems generally have demographic fluctuations due to birth and death processes. At the same time, they are exposed to changing environments. We studied populations composed of two phenotypes of bacteria and analyzed the impact that both types of fluctuations have on the mean time to extinction of the entire population if extinction is the final fate. Our results are based on Gillespie simulations and on the WKB approach applied to classical stochastic systems, here in certain limiting cases. As a function of the frequency of environmental changes, we observe a non-monotonic dependence of the mean time to extinction. Its dependencies on other system parameters are also explored. This allows the control of the mean time to extinction to be as large or as small as possible, depending on whether extinction should be avoided or is desired from the perspective of bacteria or the perspective of hosts to which the bacteria are deleterious.

Abstract High levels of allelic diversity and strong linkage disequilibrium are found in the major histocompatibility (MHC) system in humans and other vertebrates. This article proposes several descriptive statistics that quantify the extent and pattern of strong linkage disequilibrium between pairs of highly polymorphic loci. It also develops an approximate analytic theory incorporating the effects of balancing selection, mutation, recombination, and genetic drift at two closely linked loci and compares the theoretical predictions with published surveys of the MHC class II loci, DQA1 and DQB1, in humans and nonhuman primates. The descriptive statistics proposed include the fraction of complementary haplotypes (haplotypes with D″ = 1), the fraction of excess haplotypes, and the numbers of alleles at each locus in complementary haplotypes with one or more alleles at the other locus. The model assumes the infinite alleles model of mutation and the symmetric overdominance model of selection. Analytic approximations in some cases are obtained in the strong selection, weak mutation (SSWM) limit introduced by J. Gillespie. The predictions of the approximate analysis are confirmed by simulation. Both the analytic theory and simulations show that relatively few haplotypes will be found when selection is strong and recombination is weak relative to genetic drift. The model can reproduce many of the observed patterns at DQA1 and DQB1 provided that the recombination rate is assumed to be very small.

This paper is concern with modeling cholera epidemic. Despite the advances made in understanding this disease and its treatment, cholera continues to be a major public health problem in many countries. Deterministic and stochastic models emerged in modeling of cholera epidemic, in order to understand the mechanism by which cholera disease spread, conditions for cholera disease to have minor and large outbreaks. We formulate a continuous time Markov chain model for cholera epidemic transmission from the deterministic model. The basic reproduction number (R0) and the extinction thresholds of corresponding cholera continuous time Markov chain model are derived under certain assumptions. We find that, the probability of extinction (no outbreak) is 1 if R0 < 1, but less than 1 if R0 > 1. We also carry out numerical simulations using Gillespie algorithm and Runge–Kutta method to generate the sample path of cholera continuous time Markov chain model and the solution of ordinary differential equation respectively. The results show that the sample path of continuous time Markov chain model fluctuates within the solution of the ordinary differential equation.

It is challenging to predict the molecular weight distribution (MWD) for a polymer with a branched architecture, though such information will significantly benefit the design and development of branched polymers with desired properties and functions. A Monte Carlo (MC) simulation method based on the Gillespie algorithm is developed to quickly compute the MWD of branched polymers formed through step-growth polymerization, with a branched polyetherimide from two backbone monomers (4,4′-bisphenol A dianhydride and m-phenylenediamine), a chain terminator (phthalic anhydride), and a branching agent (tris[4-(4-aminophenoxy)phenyl] ethane) as an example. This polymerization involves four reactions that can be all reduced to a condensation reaction between an amine group and a carboxylic anhydride group. A comparison between the MC simulation results and the predictions of the Flory-Stockmayer theory on MWD shows that the rates of the reactions are determined by the concentrations of the functional groups on the monomers involved in each reaction. It further shows that the Flory-Stockmayer theory predicts MWD well for systems below the gel point but starts to fail for systems around or above the gel point. However, for all the systems, the MC method can be used to reliably predict MWD no matter if they are below or above the gel point. Even for a macroscopic system, a converging distribution can be quickly obtained through MC simulations on a system of only a few hundred to a few thousand monomers that have the same molar ratios as in the macroscopic system.

Abstract The stochastic simulation algorithm commonly known as Gillespie’s algorithm (originally derived for modelling well-mixed systems of chemical reactions) is now used ubiquitously in the modelling of biological processes in which stochastic effects play an important role. In well-mixed scenarios at the sub-cellular level it is often reasonable to assume that times between successive reaction/interaction events are exponentially distributed and can be appropriately modelled as a Markov process and hence simulated by the Gillespie algorithm. However, Gillespie’s algorithm is routinely applied to model biological systems for which it was never intended. In particular, processes in which cell proliferation is important (e.g. embryonic development, cancer formation) should not be simulated naively using the Gillespie algorithm since the history-dependent nature of the cell cycle breaks the Markov process. The variance in experimentally measured cell cycle times is far less than in an exponential cell cycle time distribution with the same mean. Here we suggest a method of modelling the cell cycle that restores the memoryless property to the system and is therefore consistent with simulation via the Gillespie algorithm. By breaking the cell cycle into a number of independent exponentially distributed stages, we can restore the Markov property at the same time as more accurately approximating the appropriate cell cycle time distributions. The consequences of our revised mathematical model are explored analytically as far as possible. We demonstrate the importance of employing the correct cell cycle time distribution by recapitulating the results from two models incorporating cellular proliferation (one spatial and one non-spatial) and demonstrating that changing the cell cycle time distribution makes quantitative and qualitative differences to the outcome of the models. Our adaptation will allow modellers and experimentalists alike to appropriately represent cellular proliferation—vital to the accurate modelling of many biological processes—whilst still being able to take advantage of the power and efficiency of the popular Gillespie algorithm.

Background: African swine fever (ASF) is one of the most important foreign animal diseases to the U.S. swine industry. Stakeholders in the swine production sector are on high alert as they witness the devastation of ongoing outbreaks in some of its most important trade partner countries. Efforts to improve preparedness for ASF outbreak management are proceeding in earnest and mathematical modeling is an integral part of these efforts. Aim: This study aimed to assess the impact on within-herd transmission dynamics of African swine fever (ASF) when the models used to simulate transmission assume there is homogeneous mixing of animals within a barn. Methods: Barn-level heterogeneity was explicitly captured using a stochastic, individual pig-based, heterogeneous transmission model that considers three types of infection transmission, 1) within-pen via nose-to-nose contact; 2) between-pen via nose-to-nose contact with pigs in adjacent pens; and 3) both between- and within-pen via distance independent mechanisms (e.g., via fomites). Predictions were compared between the heterogeneous and the homogeneous Gillespie models. Results: Results showed that the predicted mean number of infectious pigs at specific time points differed greatly between the homogeneous and heterogeneous models for scenarios with low levels of between pen contacts via distance independent pathways and the differences between the two model predictions were more pronounced for the slow contact rate scenario. The heterogeneous transmission model results also showed that it may take significantly longer to detect ASF, particularly in large barns when transmission predominantly occurs via nose-to-nose contact between pigs in adjacent pens. Conclusion: The findings emphasize the need for completing preliminary explorations when working with homogeneous mixing models to ascertain their suitability to predict disease outcomes.

In this study we present simulation system based on Gillespie algorithm for generating evolutionary events in the evolution scenario of microbial population. We present Gillespie simulation system adjusted to reproducing experimental data obtained in barcoding studies – experimental techniques in microbiology allowing tracing microbial populations with very high resolution. Gillespie simulation engine is constructed by defining its state vector and rules for its modifications. In order to efficiently simulate barcoded experiment by using Gillespie algorithm we provide modification - binning cells by lineages. Different bins define components of state in the Gillespie algorithm. The elaborated simulation model captures events in microbial population growth including death, division and mutations of cells. The obtained simulation results reflect population behavior, mutation wave and mutation distribution along generations. The elaborated methodology is confronted against literature data of experimental evolution of yeast tracking clones sub-generations. Simulation model was fitted to measurements in experimental data leading to good agreement.

This article presents issues related to the evolution of media and media competences, a review and analysis of selected historical, technological and educational conditions in the context of the development of digital technologies. A comparison is also made between digital, information and media competences, current development trends and future trends. The differences and requirements between qualified media users and qualified users of information technology are becoming less and less distinct. The 3 generations of Media education - 1.0, 2.0 and 3.0 were described. The main purpose of media education in the first phase of development, referred to as media education 1.0, was to develop not only critical thinking skills towards the media and media messages, but also - in a general sense - critical attitude and autonomy. Media 2.0 education can be discussed in connection with the dynamic development of the Internet and information and communication technologies, including social media, at the beginning of the 21st century. In the scТОntТПТc dТscoursО oП rОcОnt вОars, tСО concОpt oП „alРorТtСmТc culturО” Сas appОarОd, orТРТnallв defining a set of cultural artefacts that are software products, related to video games, and now describing the phenomenon in which the Big Data logic of large-scale machine learning algorithms change how culture is practiced, processed and understood (Gillespie, 2014). This stage of evolution of Media education could be identified as Media education 3.0. AI and VR and AR can accelerate teaching and learning processes through immersion, collaboration among users, realistic simulations and multi-channel communication. The topic is quite important and current in the context of changes in the education system at various levels and the challenges involved in preparing new programs.

To date, discrete event stochastic simulations of large scale biological reaction systems are extremely compute-intensive and time-consuming. Besides, it has been widely accepted that spatial factor plays a critical role in the dynamics of most biological reaction systems. The NSM (the Next Sub-Volume Method), a spatial variation of the Gillespie’s stochastic simulation algorithm (SSA), has been proposed for spatially stochastic simulation of those systems. While being able to explore high degree of parallelism in systems, NSM is inherently sequential, which still suffers from the problem of low simulation speed. Fine-grained parallel execution is an elegant way to speed up sequential simulations. Thus, based on the discrete event simulation framework JAMES II, we design and implement a PDES (Parallel Discrete Event Simulation) TW (time warp) simulator to enable the fine-grained parallel execution of spatial stochastic simulations of biological reaction systems using the ANSM (the Abstract NSM), a parallel variation of the NSM. The simulation results of classical Lotka-Volterra biological reaction system show that our time warp simulator obtains remarkable parallel speed-up against sequential execution of the NSM.I.IntroductionThe goal of Systems biology is to obtain system-level investigations of the structure and behavior of biological reaction systems by integrating biology with system theory, mathematics and computer science [1][3], since the isolated knowledge of parts can not explain the dynamics of a whole system. As the complement of “wet-lab” experiments, stochastic simulation, being called the “dry-computational” experiment, plays a more and more important role in computing systems biology [2]. Among many methods explored in systems biology, discrete event stochastic simulation is of greatly importance [4][5][6], since a great number of researches have present that stochasticity or “noise” have a crucial effect on the dynamics of small population biological reaction systems [4][7]. Furthermore, recent research shows that the stochasticity is not only important in biological reaction systems with small population but also in some moderate/large population systems [7].To date, Gillespie’s SSA [8] is widely considered to be the most accurate way to capture the dynamics of biological reaction systems instead of traditional mathematical method [5][9]. However, SSA-based stochastic simulation is confronted with two main challenges: Firstly, this type of simulation is extremely time-consuming, since when the types of species and the number of reactions in the biological system are large, SSA requires a huge amount of steps to sample these reactions; Secondly, the assumption that the systems are spatially homogeneous or well-stirred is hardly met in most real biological systems and spatial factors play a key role in the behaviors of most real biological systems [19][20][21][22][23][24]. The next sub-volume method (NSM) [18], presents us an elegant way to access the special problem via domain partition. To our disappointment, sequential stochastic simulation with the NSM is still very time-consuming, and additionally introduced diffusion among neighbor sub-volumes makes things worse. Whereas, the NSM explores a very high degree of parallelism among sub-volumes, and parallelization has been widely accepted as the most meaningful way to tackle the performance bottleneck of sequential simulations [26][27]. Thus, adapting parallel discrete event simulation (PDES) techniques to discrete event stochastic simulation would be particularly promising. Although there are a few attempts have been conducted [29][30][31], research in this filed is still in its infancy and many issues are in need of further discussion. The next section of the paper presents the background and related work in this domain. In section III, we give the details of design and implementation of model interfaces of LP paradigm and the time warp simulator based on the discrete event simulation framework JAMES II; the benchmark model and experiment results are shown in Section IV; in the last section, we conclude the paper with some future work.II. Background and Related WorkA. Parallel Discrete Event Simulation (PDES)The notion Logical Process (LP) is introduced to PDES as the abstract of the physical process [26], where a system consisting of many physical processes is usually modeled by a set of LP. LP is regarded as the smallest unit that can be executed in PDES and each LP holds a sub-partition of the whole system’s state variables as its private ones. When a LP processes an event, it can only modify the state variables of its own. If one LP needs to modify one of its neighbors’ state variables, it has to schedule an event to the target neighbor. That is to say event message exchanging is the only way that LPs interact with each other. Because of the data dependences or interactions among LPs, synchronization protocols have to be introduced to PDES to guarantee the so-called local causality constraint (LCC) [26]. By now, there are a larger number of synchronization algorithms have been proposed, e.g. the null-message [26], the time warp (TW) [32], breath time warp (BTW) [33] and etc. According to whether can events of LPs be processed optimistically, they are generally divided into two types: conservative algorithms and optimistic algorithms. However, Dematté and Mazza have theoretically pointed out the disadvantages of pure conservative parallel simulation for biochemical reaction systems [31]. B. NSM and ANSM The NSM is a spatial variation of Gillespie’ SSA, which integrates the direct method (DM) [8] with the next reaction method (NRM) [25]. The NSM presents us a pretty good way to tackle the aspect of space in biological systems by partitioning a spatially inhomogeneous system into many much more smaller “homogeneous” ones, which can be simulated by SSA separately. However, the NSM is inherently combined with the sequential semantics, and all sub-volumes share one common data structure for events or messages. Thus, directly parallelization of the NSM may be confronted with the so-called boundary problem and high costs of synchronously accessing the common data structure [29]. In order to obtain higher efficiency of parallel simulation, parallelization of NSM has to firstly free the NSM from the sequential semantics and secondly partition the shared data structure into many “parallel” ones. One of these is the abstract next sub-volume method (ANSM) [30]. In the ANSM, each sub-volume is modeled by a logical process (LP) based on the LP paradigm of PDES, where each LP held its own event queue and state variables (see Fig. 1). In addition, the so-called retraction mechanism was introduced in the ANSM too (see algorithm 1). Besides, based on the ANSM, Wang etc. [30] have experimentally tested the performance of several PDES algorithms in the platform called YH-SUPE [27]. However, their platform is designed for general simulation applications, thus it would sacrifice some performance for being not able to take into account the characteristics of biological reaction systems. Using the similar ideas of the ANSM, Dematté and Mazza have designed and realized an optimistic simulator. However, they processed events in time-stepped manner, which would lose a specific degree of precisions compared with the discrete event manner, and it is very hard to transfer a time-stepped simulation to a discrete event one. In addition, Jeschke etc.[29] have designed and implemented a dynamic time-window simulator to execution the NSM in parallel on the grid computing environment, however, they paid main attention on the analysis of communication costs and determining a better size of the time-window.Fig. 1: the variations from SSA to NSM and from NSM to ANSMC. JAMES II JAMES II is an open source discrete event simulation experiment framework developed by the University of Rostock in Germany. It focuses on high flexibility and scalability [11][13]. Based on the plug-in scheme [12], each function of JAMES II is defined as a specific plug-in type, and all plug-in types and plug-ins are declared in XML-files [13]. Combined with the factory method pattern JAMES II innovatively split up the model and simulator, which makes JAMES II is very flexible to add and reuse both of models and simulators. In addition, JAMES II supports various types of modelling formalisms, e.g. cellular automata, discrete event system specification (DEVS), SpacePi, StochasticPi and etc.[14]. Besides, a well-defined simulator selection mechanism is designed and developed in JAMES II, which can not only automatically choose the proper simulators according to the modeling formalism but also pick out a specific simulator from a serious of simulators supporting the same modeling formalism according to the user settings [15].III. The Model Interface and SimulatorAs we have mentioned in section II (part C), model and simulator are split up into two separate parts. Thus, in this section, we introduce the designation and implementation of model interface of LP paradigm and more importantly the time warp simulator.A. The Mod Interface of LP ParadigmJAMES II provides abstract model interfaces for different modeling formalism, based on which Wang etc. have designed and implemented model interface of LP paradigm[16]. However, this interface is not scalable well for parallel and distributed simulation of larger scale systems. In our implementation, we accommodate the interface to the situation of parallel and distributed situations. Firstly, the neighbor LP’s reference is replaced by its name in LP’s neighbor queue, because it is improper even dangerous that a local LP hold the references of other LPs in remote memory space. In addition, (pseudo-)random number plays a crucial role to obtain valid and meaningful results in stochastic simulations. However, it is still a very challenge work to find a good random number generator (RNG) [34]. Thus, in order to focus on our problems, we introduce one of the uniform RNGs of JAMES II to this model interface, where each LP holds a private RNG so that random number streams of different LPs can be independent stochastically. B. The Time Warp SimulatorBased on the simulator interface provided by JAMES II, we design and implement the time warp simulator, which contains the (master-)simulator, (LP-)simulator. The simulator works strictly as master/worker(s) paradigm for fine-grained parallel and distributed stochastic simulations. Communication costs are crucial to the performance of a fine-grained parallel and distributed simulation. Based on the Java remote method invocation (RMI) mechanism, P2P (peer-to-peer) communication is implemented among all (master-and LP-)simulators, where a simulator holds all the proxies of targeted ones that work on remote workers. One of the advantages of this communication approach is that PDES codes can be transferred to various hardwire environment, such as Clusters, Grids and distributed computing environment, with only a little modification; The other is that RMI mechanism is easy to realized and independent to any other non-Java libraries. Since the straggler event problem, states have to be saved to rollback events that are pre-processed optimistically. Each time being modified, the state is cloned to a queue by Java clone mechanism. Problem of this copy state saving approach is that it would cause loads of memory space. However, the problem can be made up by a condign GVT calculating mechanism. GVT reduction scheme also has a significant impact on the performance of parallel simulators, since it marks the highest time boundary of events that can be committed so that memories of fossils (processed events and states) less than GVT can be reallocated. GVT calculating is a very knotty for the notorious simultaneous reporting problem and transient messages problem. According to our problem, another GVT algorithm, called Twice Notification (TN-GVT) (see algorithm 2), is contributed to this already rich repository instead of implementing one of GVT algorithms in reference [26] and [28].This algorithm looks like the synchronous algorithm described in reference [26] (pp. 114), however, they are essentially different from each other. This algorithm has never stopped the simulators from processing events when GVT reduction, while algorithm in reference [26] blocks all simulators for GVT calculating. As for the transient message problem, it can be neglect in our implementation, because RMI based remote communication approach is synchronized, that means a simulator will not go on its processing until the remote the massage get to its destination. And because of this, the high-costs message acknowledgement, prevalent over many classical asynchronous GVT algorithms, is not needed anymore too, which should be constructive to the whole performance of the time warp simulator.IV. Benchmark Model and Experiment ResultsA. The Lotka-Volterra Predator-prey SystemIn our experiment, the spatial version of Lotka-Volterra predator-prey system is introduced as the benchmark model (see Fig. 2). We choose the system for two considerations: 1) this system is a classical experimental model that has been used in many related researches [8][30][31], so it is credible and the simulation results are comparable; 2) it is simple but helpful enough to test the issues we are interested in. The space of predator-prey System is partitioned into a2D NXNgrid, whereNdenotes the edge size of the grid. Initially the population of the Grass, Preys and Predators are set to 1000 in each single sub-volume (LP). In Fig. 2,r1,r2,r3stand for the reaction constants of the reaction 1, 2 and 3 respectively. We usedGrass,dPreyanddPredatorto stand for the diffusion rate of Grass, Prey and Predator separately. Being similar to reference [8], we also take the assumption that the population of the grass remains stable, and thusdGrassis set to zero.R1:Grass + Prey ->2Prey(1)R2:Predator +Prey -> 2Predator(2)R3:Predator -> NULL(3)r1=0.01; r2=0.01; r3=10(4)dGrass=0.0;dPrey=2.5;dPredato=5.0(5)Fig. 2: predator-prey systemB. Experiment ResultsThe simulation runs have been executed on a Linux Cluster with 40 computing nodes. Each computing node is equipped with two 64bit 2.53 GHz Intel Xeon QuadCore Processors with 24GB RAM, and nodes are interconnected with Gigabit Ethernet connection. The operating system is Kylin Server 3.5, with kernel 2.6.18. Experiments have been conducted on the benchmark model of different size of mode to investigate the execution time and speedup of the time warp simulator. As shown in Fig. 3, the execution time of simulation on single processor with 8 cores is compared. The result shows that it will take more wall clock time to simulate much larger scale systems for the same simulation time. This testifies the fact that larger scale systems will leads to more events in the same time interval. More importantly, the blue line shows that the sequential simulation performance declines very fast when the mode scale becomes large. The bottleneck of sequential simulator is due to the costs of accessing a long event queue to choose the next events. Besides, from the comparison between group 1 and group 2 in this experiment, we could also conclude that high diffusion rate increased the simulation time greatly both in sequential and parallel simulations. This is because LP paradigm has to split diffusion into two processes (diffusion (in) and diffusion (out) event) for two interactive LPs involved in diffusion and high diffusion rate will lead to high proportional of diffusion to reaction. In the second step shown in Fig. 4, the relationship between the speedups from time warp of two different model sizes and the number of work cores involved are demonstrated. The speedup is calculated against the sequential execution of the spatial reaction-diffusion systems model with the same model size and parameters using NSM.Fig. 4 shows the comparison of speedup of time warp on a64X64grid and a100X100grid. In the case of a64X64grid, under the condition that only one node is used, the lowest speedup (a little bigger than 1) is achieved when two cores involved, and the highest speedup (about 6) is achieved when 8 cores involved. The influence of the number of cores used in parallel simulation is investigated. In most cases, large number of cores could bring in considerable improvements in the performance of parallel simulation. Also, compared with the two results in Fig. 4, the simulation of larger model achieves better speedup. Combined with time tests (Fig. 3), we find that sequential simulator’s performance declines sharply when the model scale becomes very large, which makes the time warp simulator get better speed-up correspondingly.Fig. 3: Execution time (wall clock time) of Seq. and time warp with respect to different model sizes (N=32, 64, 100, and 128) and model parameters based on single computing node with 8 cores. Results of the test are grouped by the diffusion rates (Group 1: Sequential 1 and Time Warp 1. dPrey=2.5, dPredator=5.0; Group 2: dPrey=0.25, dPredator=0.5, Sequential 2 and Time Warp 2).Fig. 4: Speedup of time warp with respect to the number of work cores and the model size (N=64 and 100). Work cores are chose from one computing node. Diffusion rates are dPrey=2.5, dPredator=5.0 and dGrass=0.0.V. Conclusion and Future WorkIn this paper, a time warp simulator based on the discrete event simulation framework JAMES II is designed and implemented for fine-grained parallel and distributed discrete event spatial stochastic simulation of biological reaction systems. Several challenges have been overcome, such as state saving, roll back and especially GVT reduction in parallel execution of simulations. The Lotka-Volterra Predator-Prey system is chosen as the benchmark model to test the performance of our time warp simulator and the best experiment results show that it can obtain about 6 times of speed-up against the sequential simulation. The domain this paper concerns with is in the infancy, many interesting issues are worthy of further investigated, e.g. there are many excellent PDES optimistic synchronization algorithms (e.g. the BTW) as well. Next step, we would like to fill some of them into JAMES II. In addition, Gillespie approximation methods (tau-leap[10] etc.) sacrifice some degree of precision for higher simulation speed, but still could not address the aspect of space of biological reaction systems. The combination of spatial element and approximation methods would be very interesting and promising; however, the parallel execution of tau-leap methods should have to overcome many obstacles on the road ahead.AcknowledgmentThis work is supported by the National Natural Science Foundation of China (NSF) Grant (No.60773019) and the Ph.D. Programs Foundation of Ministry of Education of China (No. 200899980004). The authors would like to show their great gratitude to Dr. Jan Himmelspach and Dr. Roland Ewald at the University of Rostock, Germany for their invaluable advice and kindly help with JAMES II.ReferencesH. Kitano, "Computational systems biology." Nature, vol. 420, no. 6912, pp. 206-210, November 2002.H. Kitano, "Systems biology: a brief overview." Science (New York, N.Y.), vol. 295, no. 5560, pp. 1662-1664, March 2002.A. Aderem, "Systems biology: Its practice and challenges," Cell, vol. 121, no. 4, pp. 511-513, May 2005. [Online]. Available: http://dx.doi.org/10.1016/j.cell.2005.04.020.H. de Jong, "Modeling and simulation of genetic regulatory systems: A literature review," Journal of Computational Biology, vol. 9, no. 1, pp. 67-103, January 2002.C. W. Gardiner, Handbook of Stochastic Methods: for Physics, Chemistry and the Natural Sciences (Springer Series in Synergetics), 3rd ed. Springer, April 2004.D. T. Gillespie, "Simulation methods in systems biology," in Formal Methods for Computational Systems Biology, ser. Lecture Notes in Computer Science, M. Bernardo, P. Degano, and G. Zavattaro, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, vol. 5016, ch. 5, pp. 125-167.Y. Tao, Y. Jia, and G. T. Dewey, "Stochastic fluctuations in gene expression far from equilibrium: Omega expansion and linear noise approximation," The Journal of Chemical Physics, vol. 122, no. 12, 2005.D. T. Gillespie, "Exact stochastic simulation of coupled chemical reactions," Journal of Physical Chemistry, vol. 81, no. 25, pp. 2340-2361, December 1977.D. T. Gillespie, "Stochastic simulation of chemical kinetics," Annual Review of Physical Chemistry, vol. 58, no. 1, pp. 35-55, 2007.D. T. Gillespie, "Approximate accelerated stochastic simulation of chemically reacting systems," The Journal of Chemical Physics, vol. 115, no. 4, pp. 1716-1733, 2001.J. Himmelspach, R. Ewald, and A. M. Uhrmacher, "A flexible and scalable experimentation layer," in WSC '08: Proceedings of the 40th Conference on Winter Simulation. Winter Simulation Conference, 2008, pp. 827-835.J. Himmelspach and A. M. Uhrmacher, "Plug'n simulate," in 40th Annual Simulation Symposium (ANSS'07). Washington, DC, USA: IEEE, March 2007, pp. 137-143.R. Ewald, J. Himmelspach, M. Jeschke, S. Leye, and A. M. Uhrmacher, "Flexible experimentation in the modeling and simulation framework james ii-implications for computational systems biology," Brief Bioinform, vol. 11, no. 3, pp. bbp067-300, January 2010.A. Uhrmacher, J. Himmelspach, M. Jeschke, M. John, S. Leye, C. Maus, M. Röhl, and R. Ewald, "One modelling formalism & simulator is not enough! a perspective for computational biology based on james ii," in Formal Methods in Systems Biology, ser. Lecture Notes in Computer Science, J. Fisher, Ed. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, vol. 5054, ch. 9, pp. 123-138. [Online]. Available: http://dx.doi.org/10.1007/978-3-540-68413-8_9.R. Ewald, J. Himmelspach, and A. M. Uhrmacher, "An algorithm selection approach for simulation systems," pads, vol. 0, pp. 91-98, 2008.Bing Wang, Jan Himmelspach, Roland Ewald, Yiping Yao, and Adelinde M Uhrmacher. Experimental analysis of logical process simulation algorithms in james ii[C]// In M. D. Rossetti, R. R. Hill, B. Johansson, A. Dunkin, and R. G. Ingalls, editors, Proceedings of the Winter Simulation Conference, IEEE Computer Science, 2009. 1167-1179.Ewald, J. Rössel, J. Himmelspach, and A. M. Uhrmacher, "A plug-in-based architecture for random number generation in simulation systems," in WSC '08: Proceedings of the 40th Conference on Winter Simulation. Winter Simulation Conference, 2008, pp. 836-844.J. Elf and M. Ehrenberg, "Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases." Systems biology, vol. 1, no. 2, pp. 230-236, December 2004.K. Takahashi, S. Arjunan, and M. Tomita, "Space in systems biology of signaling pathways? Towards intracellular molecular crowding in silico," FEBS Letters, vol. 579, no. 8, pp. 1783-1788, March 2005.J. V. Rodriguez, J. A. Kaandorp, M. Dobrzynski, and J. G. Blom, "Spatial stochastic modelling of the phosphoenolpyruvate-dependent phosphotransferase (pts) pathway in escherichia coli," Bioinformatics, vol. 22, no. 15, pp. 1895-1901, August 2006.D. Ridgway, G. Broderick, and M. Ellison, "Accommodating space, time and randomness in network simulation," Current Opinion in Biotechnology, vol. 17, no. 5, pp. 493-498, October 2006.J. V. Rodriguez, J. A. Kaandorp, M. Dobrzynski, and J. G. Blom, "Spatial stochastic modelling of the phosphoenolpyruvate-dependent phosphotransferase (pts) pathway in escherichia coli," Bioinformatics, vol. 22, no. 15, pp. 1895-1901, August 2006.W. G. Wilson, A. M. Deroos, and E. Mccauley, "Spatial instabilities within the diffusive lotka-volterra system: Individual-based simulation results," Theoretical Population Biology, vol. 43, no. 1, pp. 91-127, February 1993.K. Kruse and J. Elf. Kinetics in spatially extended systems. In Z. Szallasi, J. Stelling, and V. Periwal, editors, System Modeling in Cellular Biology. From Concepts to Nuts and Bolts, pages 177–198. MIT Press, Cambridge, MA, 2006.M. A. Gibson and J. Bruck, "Efficient exact stochastic simulation of chemical systems with many species and many channels," The Journal of Physical Chemistry A, vol. 104, no. 9, pp. 1876-1889, March 2000.R. M. Fujimoto, Parallel and Distributed Simulation Systems (Wiley Series on Parallel and Distributed Computing). Wiley-Interscience, January 2000.Y. Yao and Y. Zhang, “Solution for analytic simulation based on parallel processing,” Journal of System Simulation, vol. 20, No.24, pp. 6617–6621, 2008.G. Chen and B. K. Szymanski, "Dsim: scaling time warp to 1,033 processors," in WSC '05: Proceedings of the 37th conference on Winter simulation. Winter Simulation Conference, 2005, pp. 346-355.M. Jeschke, A. Park, R. Ewald, R. Fujimoto, and A. M. Uhrmacher, "Parallel and distributed spatial simulation of chemical reactions," in 2008 22nd Workshop on Principles of Advanced and Distributed Simulation. Washington, DC, USA: IEEE, June 2008, pp. 51-59.B. Wang, Y. Yao, Y. Zhao, B. Hou, and S. Peng, "Experimental analysis of optimistic synchronization algorithms for parallel simulation of reaction-diffusion systems," High Performance Computational Systems Biology, International Workshop on, vol. 0, pp. 91-100, October 2009.L. Dematté and T. Mazza, "On parallel stochastic simulation of diffusive systems," in Computational Methods in Systems Biology, M. Heiner and A. M. Uhrmacher, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, vol. 5307, ch. 16, pp. 191-210.D. R. Jefferson, "Virtual time," ACM Trans. Program. Lang. Syst., vol. 7, no. 3, pp. 404-425, July 1985.J. S. Steinman, "Breathing time warp," SIGSIM Simul. Dig., vol. 23, no. 1, pp. 109-118, July 1993. [Online]. Available: http://dx.doi.org/10.1145/174134.158473 S. K. Park and K. W. Miller, "Random number generators: good ones are hard to find," Commun. ACM, vol. 31, no. 10, pp. 1192-1201, October 1988.

Bottom-up models of functionally relevant patterns of neural activity provide an explicit link between neuronal dynamics and computation. A prime example of functional activity patterns are propagating bursts of place-cell activities called hippocampal replay, which is critical for memory consolidation. The sudden and repeated occurrences of these burst states during ongoing neural activity suggest metastable neural circuit dynamics. As metastability has been attributed to noise and/or slow fatigue mechanisms, we propose a concise mesoscopic model which accounts for both. Crucially, our model is bottom-up: it is analytically derived from the dynamics of finite-size networks of Linear-Nonlinear Poisson neurons with short-term synaptic depression. As such, noise is explicitly linked to stochastic spiking and network size, and fatigue is explicitly linked to synaptic dynamics. To derive the mesoscopic model, we first consider a homogeneous spiking neural network and follow the temporal coarse-graining approach of Gillespie to obtain a “chemical Langevin equation”, which can be naturally interpreted as a stochastic neural mass model. The Langevin equation is computationally inexpensive to simulate and enables a thorough study of metastable dynamics in classical setups (population spikes and Up-Down-states dynamics) by means of phase-plane analysis. An extension of the Langevin equation for small network sizes is also presented. The stochastic neural mass model constitutes the basic component of our mesoscopic model for replay. We show that the mesoscopic model faithfully captures the statistical structure of individual replayed trajectories in microscopic simulations and in previously reported experimental data. Moreover, compared to the deterministic Romani-Tsodyks model of place-cell dynamics, it exhibits a higher level of variability regarding order, direction and timing of replayed trajectories, which seems biologically more plausible and could be functionally desirable. This variability is the product of a new dynamical regime where metastability emerges from a complex interplay between finite-size fluctuations and local fatigue.

Background. Relapse affects about 50% of AML patients who achieved remission after treatment, and the prognosis of relapsed AML is poor. Current evidence has shown that in many patients, mutations giving rise to relapse are already present at diagnosis and remain in small numbers in remission, defined as the minimal residual disease (MRD) [1]. Chemoresistant clones contributing to relapse of the disease arise from minimal residual disease (MRD) rather than resulting from newly acquired mutations during or after chemotherapy. MRD is the presence of measurable leukemic cells using non-morphologic assays. It is considered a strong predictor of relapse. The dynamics of clones comprising MRD is poorly understood and is considered influenced by a form of Darwinian selection.Methods. We propose a stochastic model based on a multitype (multi-clone) age-dependent Markov branching process to study how random events in MRD contribute to the heterogeneity in response to treatment in a cohort of six patients from The Cancer Genome Atlas database with whole genome sequencing data at two time points. Because human bone marrow cell counts are too large for direct stochastic simulation methods to be effective, we developed a hybrid numerical algorithm combining stochastic Gillespie-type and tau-leaping algorithms, and a deterministic differential equation solver, which uses much less computer time than a "straight Gillespie algorithm". Results. We developed a stochastic model of clonal evolution based on a multitype (multi-clone) age-dependent Markov branching process model of cell proliferation. In brief, we consider the critical time interval between diagnosis and initial relapse of AML that includes cytotoxic chemotherapy, chemotherapy-induced myelosuppression and decrease in leukemic cells, non-leukemic marrow recovery, and growth of the leukemic clones due to refractory or relapsed disease. Underlying our model are assumptions regarding the structure of growth, differentiation, and competition of the normal and leukemic clones. Our model reflects the stochasticity inherent when leukemic clones are near depletion after chemotherapy, which we hypothesize strongly contributes to the interpatient heterogeneity in treatment response. The parameters are estimated by fitting the expected-value model to the patient's clinical data. The available data at diagnosis includes patient's weight, percent cellularity, white blood cell count, percentage of blasts in both peripheral blood and bone marrow, and percentage of normal neutrophils in the peripheral blood. Importantly, the time to relapse and percentage of blasts in bone marrow at relapse are available. The parameters fitted to the expected-value model offer an explanation of how a leukemic clone can escape chemotherapy and promote relapse. These clones have either high proliferation rates or high self-renewal rates. As a result, there is a range of different parameter combinations that can explain their ability to succeed. On the other hand, we also study the clones that have been eradicated by the time of relapse and conclude that these clones might be eliminated either because they are not competitive and therefore surrender to other clones, or they are simply killed by chemotherapy. Also, we checked if the parameters are biologically relevant by using the model to compute the corresponding clonal growth rates for each patient. That these values fit in the clinically observed range independently found in for patients with the NPM1 mutations [1], suggests that the model is consistent with clinical data. Conclusions. Our model offers a more accurate understanding of how relapse arises and which properties allow a leukemic clone to thrive in the Darwinian competition among leukemic and normal hematopoietic clones. The model suggests a quantitative relationship between MRD and time to relapse and therefore may aid clinicians in determining when and how to implement treatment changes to postpone or prevent the time to relapse. [1]Assessment of Minimal Residual Disease in Standard-Risk AML. N Engl J Med. 2016;374:422. Relationship between MRD and time to relapse for all patients. Estimates of MRD and time to relapse are mean values from 1000 stochastic simulations. Each color corresponds to a single patient; left triangles, circles, and right triangles correspond to three parameter sets. Simulated points are fitted with a sigmoidal Hill function. Figure Disclosures No relevant conflicts of interest to declare.

In this paper we discuss a noisy mean field model for the genetic toggle switch. We show that this model approximates very well the characteristics of the system, observed using the exact Gillespie stochastic simulation algorithm. Also, we show that the system can be made exponentially stable depending on reaction parameters.

Jumps can be seen in many natural processes. Classical deterministic modeling approach explains the dynamical behavior of such systems by using impulsive differential equations. This modeling strategy assumes that the dynamical behavior of the whole system is deterministic, continuous, and it adds jumps to the state vector at certain times. Although deterministic approach is satisfactory in many cases, it is a well-known fact that stochasticity or uncertainty has crucial importance for dynamical behavior of many others. In this study, we propose to include this abrupt change in the stochastic modeling approach, beside the deterministic one. In our model, we introduce jumps to chemical master equation and use the Gillespie direct method to simulate the evolutionary system. To illustrate the idea and distinguish the differences, we present some numerically solved examples.

Abstract Summary DIFFpop is an R package designed to simulate cellular differentiation hierarchies using either exponentially-expanding or fixed population sizes. The software includes functionalities to simulate clonal evolution due to the emergence of driver mutations under the infinite-allele assumption as well as options for simulation and analysis of single cell barcoding and labeling data. The software uses the Gillespie Stochastic Simulation Algorithm and a modification of expanding or fixed-size stochastic process models expanded to a large number of cell types and scenarios. Availability and implementation DIFFpop is available as an R-package along with vignettes on Github (https://github.com/ferlicjl/diffpop). Supplementary information Supplementary data are available at Bioinformatics online.

Summary In this paper we study the intrinsic noise effect on the switching behavior of a simple genetic circuit corresponding to the genetic toggle switch model. The numerical results obtained from a noisy mean-field model are compared to those obtained from the stochastic Gillespie simulation of the corresponding system of chemical reactions. Our results show that by using a two step reaction approach for modeling the transcription and translation processes one can make the system to lock in one of the steady states for exponentially long times.

Active intracellular transport is mainly performed by a group of special nanomachines called motor proteins. During transport, cooperation between motor proteins significantly influences important transport features, such as distance and velocity. To understand this mechanism, we combine Gillespie simulation and analytical derivation to demonstrate how the mechanical properties of a single motor influence the cooperation between multiple motors, further regulating the transport distance. In addition, we build a deep learning model to help us quickly obtain the motor parameters. Our results shed light on the physical nature of intracellular transport by motor proteins with the same directionality.

Abstract We describe a hybrid simulation technique that uses the Nernst-Planck (NP) transport equation to compute steady-state ionic flux in a non-equilibrium system and uses the Local Equilibrium Monte Carlo (LEMC) simulation technique to establish the statistical mechanical relation between the two crucial functions present in the NP equation: the concentration and the electrochemical potential profiles (Boda, D., Gillespie, D., J. Chem. Theor. Comput., 2012 8(3), 824–829). The LEMC method is an adaptation of the Grand Canonical Monte Carlo method to a non-equilibrium situation. We apply the resulting NP+LEMC method to ionic systems, where two reservoirs of electrolytes are separated by a membrane that allows the diffusion of ions through a nanopore. The nanopore can be natural (as the calcium selective Ryanodine Receptor ion channel) or synthetic (as a rectifying bipolar nanopore). We show results for these two systems and demonstrate the effectiveness of the NP+LEMC technique.

Sequence simulators are fundamental tools in bioinformatics, as they allow us to test data processing and inference tools, and are an essential component of some inference methods. The ongoing surge in available sequence data is however testing the limits of our bioinformatics software. One example is the large number of SARS-CoV-2 genomes available, which are beyond the processing power of many methods, and simulating such large datasets is also proving difficult. Here, we present a new algorithm and software for efficiently simulating sequence evolution along extremely large trees (e.g. > 100, 000 tips) when the branches of the tree are short, as is typical in genomic epidemiology. Our algorithm is based on the Gillespie approach, and it implements an efficient multi-layered search tree structure that provides high computational efficiency by taking advantage of the fact that only a small proportion of the genome is likely to mutate at each branch of the considered phylogeny. Our open source software allows easy integration with other Python packages as well as a variety of evolutionary models, including indel models and new hypermutability models that we developed to more realistically represent SARS-CoV-2 genome evolution.

Many methods have been used to model epidemic spreading. They include ordinary differential equation systems for globally homogeneous environments and partial differential equation systems to take into account spatial localisation and inhomogeneity. Stochastic differential equations systems have been used to model the inherent stochasticity of epidemic spreading processes. In our case study, we wanted to model the numbers of individuals in different states of the disease, and their locations in the country. Among the many existing methods we used our own variant of the well known Gillespie stochastic algorithm, along with the sub-volumes method to take into account the spatial localisation. Our algorithm allows us to easily switch from stochastic discrete simulation to continuous deterministic resolution using mean values. We applied our approaches on the study of the Covid-19 epidemic in France. The stochastic discrete version of Pandæsim showed very good correlations between the simulation results and the statistics gathered from hospitals, both on day by day and on global numbers, including the effects of the lockdown. Moreover, we have highlighted interesting differences in behaviour between the continuous and discrete methods that may arise in some particular conditions.

In this paper, we present TauSSA, a discrete-event simulation tool for stochastic queueing networks integrated in the LINE solver. TauSSA combines Gillespie's stochastic simulation algorithm with tau leaping, a methodology for optimistic simulation acceleration. Although tau leaping is frequently used in chemical reaction network simulation, it has so far found limited application in queueing theory. TauSSA offers one of the very first attempts to make this method broadly applicable to analyze extended queueing network models, which include class switching, fork-join, and non-exponential service and arrival distributions. We conceptualize various strategies for handling ordering and illegal states in tau leaping that arise specifically within queueing network models, and compare their performance through numerical experiments. Our main finding is that strategies that sort events based on the network topological order incur a better trade-off between speedup and approximation error.

ONTARE. REVISTA DE INVESTIGACIÓN DE LA FACULTAD DE INGENIERÍAEste artículo se propone presentar dos de los principales enfoques que están disponibles hoy en día para modelar y simular la evolución dinámica de sistemas vivos: el modelo determinístico continuo, que es dictado por los sistemas de ecuaciones diferenciales ordinarias, y el estocástico discreto, que encuentra su base en el algoritmo de simulación estocástica propuesto por Gillespie en 1976. El objetivo de esta comparación es proporcionar la información necesaria para apoyar la selección de un enfoque de modelaje, basado en un conjunto de criterios verificables. Para alcanzar este objetivo, se analizan los fundamentos de la modelación, se propone un ejemplo de modelado para un sistema de vida simple y se discuten las principales ventajas y desventajas de cada enfoque ABSTRACT This research paper describes two main approaches that are currently available for modelling and simulating the evolution of dynamic living systems, namely the continuous-deterministic one, which is rendered by systems of ordinary differential equations and the discrete-stochastic one, which bases on the Stochastic simulation algorithm proposed by Gil/espíe in 1976. The aim of this comparison is to provide the necessary information to support the selection of a modelling approach, focusing on a set of verifiable criteria. For this reason, we review the rationale approach of modelling, proposing a modelling sample for a single living system and discussing the main advantages and drawbacks of each approach.