Stochastic simulation in systems biology

Dynamics of two feed forward genetic motifs in the presence of molecular noise

Understanding the function of common motifs in gene regulatory networks is an important goal of systems biology. Feed forward loops (FFLs) are an example of such a motif. In FFLs, a gene (X) regulates another gene (Z) both directly and via an intermediary gene (Y). Previous theoretical studies have suggested several possible functions for FFLs, based on their transient responses to changes in input signals (using deterministic models) and their fluctuations around steady state (using stochastic models). In this paper we study stochastic models of the two most common FFLs, "coherent type 1" and "incoherent type 1". We incorporate molecular noise by treating DNA binding, transcription, translation, and decay as stochastic processes. By comparing the dynamics of these loops with models of alternative networks (in which X does not regulate Y), we explore how FFLs act to process information in the presence of noise. This work highlights the importance of incorporating realistic molecular noise in studying both the transient and steady-state behavior of gene regulatory networks.

Application of a novel numerical simulation to biochemical reaction systems

Recent advancements in omics and single-cell analysis highlight the necessity of numerical methods for managing the complexity of biological data. This paper introduces a simulation program for biochemical reaction systems based on the natural number simulation (NNS) method. This novel approach ensures the equitable treatment of all molecular entities, such as DNA, proteins, H2O, and hydrogen ions (H+), in biological systems. Central to NNS is its use of stoichiometric formulas, simplifying the modeling process and facilitating efficient and accurate simulations of diverse biochemical reactions. The advantage of this method is its ability to manage all molecules uniformly, ensuring a balanced representation in simulations. Detailed in Python, NNS is adept at simulating various reactions, ranging from water ionization to Michaelis-Menten kinetics and complex gene-based systems, making it an effective tool for scientific and engineering research.

Holimap: an accurate and efficient method for solving stochastic gene network dynamics

Gene-gene interactions are crucial to the control of sub-cellular processes but our understanding of their stochastic dynamics is hindered by the lack of simulation methods that can accurately and efficiently predict how the distributions of gene product numbers vary across parameter space. To overcome these difficulties, here we present Holimap (high-order linear-mapping approximation), an approach that approximates the protein or mRNA number distributions of a complex gene regulatory network by the distributions of a much simpler reaction system. We demonstrate Holimap's computational advantages over conventional methods by applying it to predict the stochastic time-dependent dynamics of various gene networks, including transcriptional networks ranging from simple autoregulatory loops to complex randomly connected networks, post-transcriptional networks, and post-translational networks. Holimap is ideally suited to study how the intricate network of gene-gene interactions results in precise coordination and control of gene expression.

Revisiting kinetic Monte Carlo algorithms for time-dependent processes: From open-loop control to feedback control

Simulating stochastic systems with feedback control is challenging due to the complex interplay between the system's dynamics and the feedback-dependent control protocols. We present a single-step-trajectory probability analysis to time-dependent stochastic systems. Based on this analysis, we revisit several time-dependent kinetic Monte Carlo (KMC) algorithms designed for systems under open-loop-control protocols. Our analysis provides a unified alternative proof to these algorithms, summarized into a pedagogical tutorial. Moreover, with the trajectory probability analysis, we present a novel feedback-controlled KMC algorithm that accurately captures the dynamics systems controlled by an external signal based on the measurements of the system's state. Our method correctly captures the system dynamics and avoids the artificial Zeno effect that arises from incorrectly applying the direct Gillespie algorithm to feedback-controlled systems. This work provides a unified perspective on existing open-loop-control KMC algorithms and also offers a powerful and accurate tool for simulating stochastic systems with feedback control.

Virtual patient with temporal evolution for mechanical ventilation trial studies: A stochastic model approach

BACKGROUND AND OBJECTIVE

Healthcare datasets are plagued by issues of data scarcity and class imbalance. Clinically validated virtual patient (VP) models can provide accurate in-silico representations of real patients and thus a means for synthetic data generation in hospital critical care settings. This research presents a realistic, time-varying mechanically ventilated respiratory failure VP profile synthesised using a stochastic model.

METHODS

A stochastic model was developed using respiratory elastance (Ers) data from two clinical cohorts and averaged over 30-minute time intervals. The stochastic model was used to generate future Ers data based on current Ers values with added normally distributed random noise. Self-validation of the VPs was performed via Monte Carlo simulation and retrospective Ers profile fitting. A stochastic VP cohort of temporal Ers evolution was synthesised and then compared to an independent retrospective patient cohort data in a virtual trial across several measured patient responses, where similarity of profiles validates the realism of stochastic model generated VP profiles.

RESULTS

A total of 120,000 3-hour VPs for pressure control (PC) and volume control (VC) ventilation modes are generated using stochastic simulation. Optimisation of the stochastic simulation process yields an ideal noise percentage of 5-10% and simulation iteration of 200,000 iterations, allowing the simulation of a realistic and diverse set of Ers profiles. Results of self-validation show the retrospective Ers profiles were able to be recreated accurately with a mean squared error of only 0.099 [0.009-0.790]% for the PC cohort and 0.051 [0.030-0.126]% for the VC cohort. A virtual trial demonstrates the ability of the stochastic VP cohort to capture Ers trends within and beyond the retrospective patient cohort providing cohort-level validation.

CONCLUSION

VPs capable of temporal evolution demonstrate feasibility for use in designing, developing, and optimising bedside MV guidance protocols through in-silico simulation and validation. Overall, the temporal VPs developed using stochastic simulation alleviate the need for lengthy, resource intensive, high cost clinical trials, while facilitating statistically robust virtual trials, ultimately leading to improved patient care and outcomes in mechanical ventilation.

Fluorescent biosensor imaging meets deterministic mathematical modelling: quantitative investigation of signalling compartmentalization

Cells execute specific responses to diverse environmental cues by encoding information in distinctly compartmentalized biochemical signalling reactions. Genetically encoded fluorescent biosensors enable the spatial and temporal monitoring of signalling events in live cells. Temporal and spatiotemporal computational models can be used to interpret biosensor experiments in complex biochemical networks and to explore hypotheses that are difficult to test experimentally. In this review, we first provide brief discussions of the experimental toolkit of fluorescent biosensors as well as computational basics with a focus on temporal and spatiotemporal deterministic models. We then describe how we used this combined approach to identify and investigate a protein kinase A (PKA) - cAMP - Ca2+ oscillatory circuit in MIN6 β cells, a mouse pancreatic β cell system. We describe the application of this combined approach to interrogate how this oscillatory circuit is differentially regulated in a nano-compartment formed at the plasma membrane by the scaffolding protein A kinase anchoring protein 79/150. We leveraged both temporal and spatiotemporal deterministic models to identify the key regulators of this oscillatory circuit, which we confirmed with further experiments. The powerful approach of combining live-cell biosensor imaging with quantitative modelling, as discussed here, should find widespread use in the investigation of spatiotemporal regulation of cell signalling.

Predicting mechanically ventilated patients future respiratory system elastance - A stochastic modelling approach

BACKGROUND AND OBJECTIVE

Respiratory mechanics of mechanically ventilated patients evolve significantly with time, disease state and mechanical ventilation (MV) treatment. Existing deterministic data prediction methods fail to comprehensively describe the multiple sources of heterogeneity of biological systems. This research presents two respiratory mechanics stochastic models with increased prediction accuracy and range, offering improved clinical utility in MV treatment.

METHODS

Two stochastic models (SM2 and SM3) were developed using retrospective patient respiratory elastance (E_rs) from two clinical cohorts which were averaged over time intervals of 10 and 30 min respectively. A stochastic model from a previous study (SM1) was used to benchmark performance. The stochastic models were clinically validated on an independent retrospective clinical cohort of 14 patients. Differences in predictive ability were evaluated using the difference in percentile lines and cumulative distribution density (CDD) curves.

RESULTS

Clinical validation shows all three models captured more than 98% (median) of future E_rs data within the 5th - 95th percentile range. Comparisons of stochastic model percentile lines reported a maximum mean absolute percentage difference of 5.2%. The absolute differences of CDD curves were less than 0.25 in the ranges of 5 < E_rs (cmH₂O/L) < 85, suggesting similar predictive capabilities within this clinically relevant E_rs range.

CONCLUSION

The new stochastic models significantly improve prediction, clinical utility, and thus feasibility for synchronisation with clinical interventions. Paired with other MV protocols, the stochastic models developed can potentially form part of decision support systems, providing guided, personalised, and safe MV treatment.

A computational design of a programmable biological processor

Basic synthetic information processing structures, such as logic gates, oscillators and flip-flops, have already been implemented in living organisms. Current implementations of these structures have yet to be extended to more complex processing structures that would constitute a biological computer. We make a step forward towards the construction of a biological computer. We describe a model-based computational design of a biological processor that uses transcription and translation resources of the host cell to perform its operations. The proposed processor is composed of an instruction memory containing a biological program, a program counter that is used to address this memory, and a biological oscillator that triggers the execution of the next instruction in the memory. We additionally describe the implementation of a biological compiler that compiles a sequence of human-readable instructions into ordinary differential equation-based models, which can be used to simulate and analyse the dynamics of the processor. The proposed implementation presents the first programmable biological processor that exploits cellular resources to execute the specified instructions. We demonstrate the application of the described processor on a set of simple yet scalable biological programs. Biological descriptions of these programs can be produced manually or automatically using the provided compiler.

Bayesian Inference of Stochastic Dynamic Models Using Early-Rejection Methods Based on Sequential Stochastic Simulations

Stochastic modelling is an important method to investigate the functions of noise in a wide range of biological systems. However, the parameter inference for stochastic models is still a challenging problem partially due to the large computing time required for stochastic simulations. To address this issue, we propose a novel early-rejection method by using sequential stochastic simulations. We first show that a large number of stochastic simulations are required to obtain reliable inference results. Instead of generating a large number of simulations for each parameter sample, we propose to generate these simulations in a number of stages. The simulation process will go to the next stage only if the accuracy of simulations at the current stage satisfies a given error criterion. We propose a formula to determine the error criterion and use a stochastic differential equation model to examine the effects of different criteria. Three biochemical network models are used to evaluate the efficiency and accuracy of the proposed method. Numerical results suggest the proposed early-rejection method achieves substantial improvement in the efficiency for the inference of stochastic models.

Multiscale Computational Modeling of Vascular Adaptation: A Systems Biology Approach Using Agent-Based Models

The widespread incidence of cardiovascular diseases and associated mortality and morbidity, along with the advent of powerful computational resources, have fostered an extensive research in computational modeling of vascular pathophysiology field and promoted in-silico models as a support for biomedical research. Given the multiscale nature of biological systems, the integration of phenomena at different spatial and temporal scales has emerged to be essential in capturing mechanobiological mechanisms underlying vascular adaptation processes. In this regard, agent-based models have demonstrated to successfully embed the systems biology principles and capture the emergent behavior of cellular systems under different pathophysiological conditions. Furthermore, through their modular structure, agent-based models are suitable to be integrated with continuum-based models within a multiscale framework that can link the molecular pathways to the cell and tissue levels. This can allow improving existing therapies and/or developing new therapeutic strategies. The present review examines the multiscale computational frameworks of vascular adaptation with an emphasis on the integration of agent-based approaches with continuum models to describe vascular pathophysiology in a systems biology perspective. The state-of-the-art highlights the current gaps and limitations in the field, thus shedding light on new areas to be explored that may become the future research focus. The inclusion of molecular intracellular pathways (e.g., genomics or proteomics) within the multiscale agent-based modeling frameworks will certainly provide a great contribution to the promising personalized medicine. Efforts will be also needed to address the challenges encountered for the verification, uncertainty quantification, calibration and validation of these multiscale frameworks.

A mathematical model which examines age-related stochastic fluctuations in DNA maintenance methylation

Due to its complexity and its ubiquitous nature the ageing process remains an enduring biological puzzle. Many molecular mechanisms and biochemical process have become synonymous with ageing. However, recent findings have pinpointed epigenetics as having a key role in ageing and healthspan. In particular age related changes to DNA methylation offer the possibility of monitoring the trajectory of biological ageing and could even be used to predict the onset of diseases such as cancer, Alzheimer's disease and cardiovascular disease. At the molecular level emerging evidence strongly suggests the regulatory processes which govern DNA methylation are subject to intracellular stochasticity. It is challenging to fully understand the impact of stochasticity on DNA methylation levels at the molecular level experimentally. An ideal solution is to use mathematical models to capture the essence of the stochasticity and its outcomes. In this paper we present a novel stochastic model which accounts for specific methylation levels within a gene promoter. Uncertainty of the eventual site-specific methylation levels for different values of methylation age, depending on the initial methylation levels were analysed. Our model predicts the observed bistable levels in CpG islands. In addition, simulations with various levels of noise indicate that uncertainty predominantly spreads through the hypermethylated region of stability, especially for large values of input noise. A key outcome of the model is that CpG islands with high to intermediate methylation levels tend to be more susceptible to dramatic DNA methylation changes due to increasing methylation age.

A Petri nets-based framework for whole-cell modeling

Whole-cell modeling aims to incorporate all main genes and processes, and their interactions of a cell in one model. Whole-cell modeling has been regarded as the central aim of systems biology but also as a grand challenge, which plays essential roles in current and future systems biology. In this paper, we analyze whole-cell modeling challenges and requirements and classify them into three aspects (or dimensions): heterogeneous biochemical networks, uncertainties in components, and representation of cell structure. We then explore how to use different Petri net classes to address different aspects of whole-cell modeling requirements. Based on these analyses, we present a Petri nets-based framework for whole-cell modeling, which not only addresses many whole-cell modeling requirements, but also offers a graphical, modular, and hierarchical modeling tool. We think this framework can offer a feasible modeling approach for whole-cell model construction.

FiCoS: A fine-grained and coarse-grained GPU-powered deterministic simulator for biochemical networks

Mathematical models of biochemical networks can largely facilitate the comprehension of the mechanisms at the basis of cellular processes, as well as the formulation of hypotheses that can be tested by means of targeted laboratory experiments. However, two issues might hamper the achievement of fruitful outcomes. On the one hand, detailed mechanistic models can involve hundreds or thousands of molecular species and their intermediate complexes, as well as hundreds or thousands of chemical reactions, a situation generally occurring in rule-based modeling. On the other hand, the computational analysis of a model typically requires the execution of a large number of simulations for its calibration, or to test the effect of perturbations. As a consequence, the computational capabilities of modern Central Processing Units can be easily overtaken, possibly making the modeling of biochemical networks a worthless or ineffective effort. To the aim of overcoming the limitations of the current state-of-the-art simulation approaches, we present in this paper FiCoS, a novel "black-box" deterministic simulator that effectively realizes both a fine-grained and a coarse-grained parallelization on Graphics Processing Units. In particular, FiCoS exploits two different integration methods, namely, the Dormand-Prince and the Radau IIA, to efficiently solve both non-stiff and stiff systems of coupled Ordinary Differential Equations. We tested the performance of FiCoS against different deterministic simulators, by considering models of increasing size and by running analyses with increasing computational demands. FiCoS was able to dramatically speedup the computations up to 855×, showing to be a promising solution for the simulation and analysis of large-scale models of complex biological processes.

A mathematical modelling framework for the regulation of intra-cellular OCT4 in human pluripotent stem cells

Human pluripotent stem cells (hPSCs) have the potential to differentiate into all cell types, a property known as pluripotency. A deeper understanding of how pluripotency is regulated is required to assist in controlling pluripotency and differentiation trajectories experimentally. Mathematical modelling provides a non-invasive tool through which to explore, characterise and replicate the regulation of pluripotency and the consequences on cell fate. Here we use experimental data of the expression of the pluripotency transcription factor OCT4 in a growing hPSC colony to develop and evaluate mathematical models for temporal pluripotency regulation. We consider fractional Brownian motion and the stochastic logistic equation and explore the effects of both additive and multiplicative noise. We illustrate the use of time-dependent carrying capacities and the introduction of Allee effects to the stochastic logistic equation to describe cell differentiation. We conclude both methods adequately capture the decline in OCT4 upon differentiation, but the Allee effect model has the advantage of allowing differentiation to occur stochastically in a sub-set of cells. This mathematical framework for describing intra-cellular OCT4 regulation can be extended to other transcription factors and developed into predictive models.

Deciphering the Role of Fluctuation Dependent Intercellular Communication in Neural Stem Cell Development

Neural stem cells (NPCs) efficiently communicate in an intercellular manner to govern specific cell fate decisions during the developmental process despite withstanding the fluctuating cellular environment. How these fluctuations from diverse origins functionally affect the precise cell fate decision making remains elusive. By taking a stochastic mathematical modeling approach, we unravel that the transcriptional variability arising within an NPC population due to intermittent cell cycle events significantly influences the neuron to NPC ratio during development. Our model proficiently quantifies the impact of different sources of heterogeneities in maintaining an exact neuron to NPC ratio and predicts plausible experimental ways to fine-tune the development of NPCs. In the future, these modeling insights may lead to better therapeutic avenues to regenerate neurons from NPCs.

Embracing Noise in Chemical Reaction Networks

We provide a short review of stochastic modeling in chemical reaction networks for mathematical and quantitative biologists. We use as case studies two publications appearing in this issue of the Bulletin, on the modeling of quasi-steady-state approximations and cell polarity. Reasons for the relevance of stochastic modeling are described along with some common differences between stochastic and deterministic models.

Towards understanding the messengers of extracellular space: Computational models of outside-in integrin reaction networks

The interactions between cells and their extracellular matrix (ECM) are critically important for homeostatic control of cell growth, proliferation, differentiation and apoptosis. Transmembrane integrin molecules facilitate the communication between ECM and the cell. Since the characterization of integrins in the late 1980s, there has been great advancement in understanding the function of integrins at different subcellular levels. However, the versatility in molecular pathways integrins are involved in, the high diversity in their interaction partners both outside and inside the cell as well as on the cell membrane and the short lifetime of events happening at the cell-ECM interface make it difficult to elucidate all the details regarding integrin function experimentally. To overcome the experimental challenges and advance the understanding of integrin biology, computational modeling tools have been used extensively. In this review, we summarize the computational models of integrin signaling while we explain the function of integrins at three main subcellular levels (outside the cell, cell membrane, cytosol). We also discuss how these computational modeling efforts can be helpful in other disciplines such as biomaterial design. As such, this review is a didactic modeling summary for biomaterial researchers interested in complementing their experimental work with computational tools or for seasoned computational scientists that would like to advance current in silico integrin models.

Identifiability analysis for stochastic differential equation models in systems biology

Mathematical models are routinely calibrated to experimental data, with goals ranging from building predictive models to quantifying parameters that cannot be measured. Whether or not reliable parameter estimates are obtainable from the available data can easily be overlooked. Such issues of parameter identifiability have important ramifications for both the predictive power of a model, and the mechanistic insight that can be obtained. Identifiability analysis is well-established for deterministic, ordinary differential equation (ODE) models, but there are no commonly adopted methods for analysing identifiability in stochastic models. We provide an accessible introduction to identifiability analysis and demonstrate how existing ideas for analysis of ODE models can be applied to stochastic differential equation (SDE) models through four practical case studies. To assess structural identifiability, we study ODEs that describe the statistical moments of the stochastic process using open-source software tools. Using practically motivated synthetic data and Markov chain Monte Carlo methods, we assess parameter identifiability in the context of available data. Our analysis shows that SDE models can often extract more information about parameters than deterministic descriptions. All code used to perform the analysis is available on Github.

Applications of personalised signalling network models in precision oncology

As our ability to provide in-depth, patient-specific characterisation of the molecular alterations within tumours rapidly improves, it is becoming apparent that new approaches will be required to leverage the power of this data and derive the full benefit for each individual patient. Systems biology approaches are beginning to emerge within this field as a potential method of incorporating large volumes of network level data and distilling a coherent, clinically-relevant prediction of drug response. However, the initial promise of this developing field is yet to be realised. Here we argue that in order to develop these precise models of individual drug response and tailor treatment accordingly, we will need to develop mathematical models capable of capturing both the dynamic nature of drug-response signalling networks and key patient-specific information such as mutation status or expression profiles. We also review the modelling approaches commonly utilised within this field, and outline recent examples of their use in furthering the application of systems biology for a precision medicine approach to cancer treatment.

Stochastic modeling of human papillomavirusearly promoter gene regulation

High risk forms of human papillomaviruses (HPVs) promote cancerous lesions and are implicated in almost all cervical cancer. Of particular relevance to cancer progression is regulation of the early promoter that controls gene expression in the initial phases of infection and can eventually lead to pre-cancer progression. Our goal was to develop a stochastic model to investigate the control mechanisms that regulate gene expression from the HPV early promoter. Our model integrates modules that account for transcriptional, post-transcriptional, translational and post-translational regulation of E1 and E2 early genes to form a functioning gene regulatory network. Each module consists of a set of biochemical steps whose stochastic evolution is governed by a chemical Master Equation and can be simulated using the Gillespie algorithm. To investigate the role of noise in gene expression, we compared our stochastic simulations with solutions to ordinary differential equations for the mean behavior of the system that are valid under the conditions of large molecular abundances and quasi-equilibrium for fast reactions. The model produced results consistent with known HPV biology. Our simulation results suggest that stochasticity plays a pivotal role in determining the dynamics of HPV gene expression. In particular, the combination of positive and negative feedback regulation generates stochastic bursts of gene expression. Analysis of the model reveals that regulation at the promoter affects burst amplitude and frequency, whereas splicing is more specialized to regulate burst frequency. Our results also suggest that splicing enhancers are a significant source of stochasticity in pre-mRNA abundance and that the number of viruses infecting the host cell represents a third important source of stochasticity in gene expression.

Estimating growth patterns and driver effects in tumor evolution from individual samples

Tumors accumulate thousands of mutations, and sequencing them has given rise to methods for finding cancer drivers via mutational recurrence. However, these methods require large cohorts and underperform for low recurrence. Recently, ultra-deep sequencing has enabled accurate measurement of VAFs (variant-allele frequencies) for mutations, allowing the determination of evolutionary trajectories. Here, based solely on the VAF spectrum for an individual sample, we report on a method that identifies drivers and quantifies tumor growth. Drivers introduce perturbations into the spectrum, and our method uses the frequency of hitchhiking mutations preceding a driver to measure this. As validation, we use simulation models and 993 tumors from the Pan-Cancer Analysis of Whole Genomes (PCAWG) Consortium with previously identified drivers. Then we apply our method to an ultra-deep sequenced acute myeloid leukemia (AML) tumor and identify known cancer genes and additional driver candidates. In summary, our framework presents opportunities for personalized driver diagnosis using sequencing data from a single individual.

Discrete and continuous models of probability flux of switching dynamics: Uncovering stochastic oscillations in a toggle-switch system

The probability flux and velocity in stochastic reaction networks can help in characterizing dynamic changes in probability landscapes of these networks. Here, we study the behavior of three different models of probability flux, namely, the discrete flux model, the Fokker-Planck model, and a new continuum model of the Liouville flux. We compare these fluxes that are formulated based on, respectively, the chemical master equation, the stochastic differential equation, and the ordinary differential equation. We examine similarities and differences among these models at the nonequilibrium steady state for the toggle switch network under different binding and unbinding conditions. Our results show that at a strong stochastic condition of weak promoter binding, continuum models of Fokker-Planck and Liouville fluxes deviate significantly from the discrete flux model. Furthermore, we report the discovery of stochastic oscillation in the toggle-switch system occurring at weak binding conditions, a phenomenon captured only by the discrete flux model.

Spatial Stochastic Intracellular Kinetics: A Review of Modelling Approaches

Models of chemical kinetics that incorporate both stochasticity and diffusion are an increasingly common tool for studying biology. The variety of competing models is vast, but two stand out by virtue of their popularity: the reaction-diffusion master equation and Brownian dynamics. In this review, we critically address a number of open questions surrounding these models: How can they be justified physically? How do they relate to each other? How do they fit into the wider landscape of chemical models, ranging from the rate equations to molecular dynamics? This review assumes no prior knowledge of modelling chemical kinetics and should be accessible to a wide range of readers.

Systems biology primer: the basic methods and approaches

Systems biology is an integrative discipline connecting the molecular components within a single biological scale and also among different scales (e.g. cells, tissues and organ systems) to physiological functions and organismal phenotypes through quantitative reasoning, computational models and high-throughput experimental technologies. Systems biology uses a wide range of quantitative experimental and computational methodologies to decode information flow from genes, proteins and other subcellular components of signaling, regulatory and functional pathways to control cell, tissue, organ and organismal level functions. The computational methods used in systems biology provide systems-level insights to understand interactions and dynamics at various scales, within cells, tissues, organs and organisms. In recent years, the systems biology framework has enabled research in quantitative and systems pharmacology and precision medicine for complex diseases. Here, we present a brief overview of current experimental and computational methods used in systems biology.

Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach

We develop numerical methods for stochastic reaction-diffusion systems based on approaches used for fluctuatinghydrodynamics (FHD). For hydrodynamicsystems, the FHD formulation is formally described by stochastic partial differential equations (SPDEs). In the reaction-diffusion systems we consider, our model becomes similar to the reaction-diffusion master equation (RDME) description when our SPDEs are spatially discretized and reactions are modeled as a source term having Poissonfluctuations. However, unlike the RDME, which becomes prohibitively expensive for an increasing number of molecules, our FHD-based description naturally extends from the regime where fluctuations are strong, i.e., each mesoscopic cell has few (reactive) molecules, to regimes with moderate or weak fluctuations, and ultimately to the deterministic limit. By treating diffusion implicitly, we avoid the severe restriction on time step size that limits all methods based on explicit treatments of diffusion and construct numerical methods that are more efficient than RDME methods, without compromising accuracy. Guided by an analysis of the accuracy of the distribution of steady-state fluctuations for the linearized reaction-diffusion model, we construct several two-stage (predictor-corrector) schemes, where diffusion is treated using a stochastic Crank-Nicolson method, and reactions are handled by the stochastic simulation algorithm of Gillespie or a weakly second-order tau leaping method. We find that an implicit midpoint tau leaping scheme attains second-order weak accuracy in the linearized setting and gives an accurate and stable structure factor for a time step size of an order of magnitude larger than the hopping time scale of diffusing molecules. We study the numerical accuracy of our methods for the Schlögl reaction-diffusion model both in and out of thermodynamic equilibrium. We demonstrate and quantify the importance of thermodynamicfluctuations to the formation of a two-dimensional Turing-like pattern and examine the effect of fluctuations on three-dimensional chemical front propagation. By comparing stochastic simulations to deterministic reaction-diffusion simulations, we show that fluctuations accelerate pattern formation in spatially homogeneous systems and lead to a qualitatively different disordered pattern behind a traveling wave.

Hybrid stochastic framework predicts efficacy of prophylaxis against HIV: An example with different dolutegravir prophylaxis schemes

To achieve the 90-90-90 goals set by UNAIDS, the number of new HIV infections needs to decrease to approximately 500,000 by 2020. One of the 'five pillars' to achieve this goal is pre-exposure prophylaxis (PrEP). Truvada (emtricitabine-tenofovir) is currently the only medication approved for PrEP. Despite its advantages, Truvada is costly and requires individuals to adhere to the once-daily regimen. To improve PrEP, many next-generation regimen, including long-acting formulations, are currently investigated. However, pre-clinical testing may not guide candidate selection, since it often fails to translate into clinical efficacy. On the other hand, quantifying prophylactic efficacy in the clinic is ethically problematic and requires to conduct long (years) and large (N>1000 individuals) trials, precluding systematic evaluation of candidates and deployment strategies. To prioritize- and help design PrEP regimen, tools are urgently needed that integrate pharmacological-, viral- and host factors determining prophylactic efficacy. Integrating the aforementioned factors, we developed an efficient and exact stochastic simulation approach to predict prophylactic efficacy, as an example for dolutegravir (DTG). Combining the population pharmacokinetics of DTG with the stochastic framework, we predicted that plasma concentrations of 145.18 and 722.23nM prevent 50- and 90% sexual transmissions respectively. We then predicted the reduction in HIV infection when DTG was used in PrEP, PrEP 'on demand' and post-exposure prophylaxis (PEP) before/after virus exposure. Once daily PrEP with 50mg oral DTG prevented 99-100% infections, and 85% of infections when 50% of dosing events were missed. PrEP 'on demand' prevented 79-84% infections and PEP >80% when initiated within 6 hours after virus exposure and continued for as long as possible. While the simulation framework can easily be adapted to other PrEP candidates, our simulations indicated that oral 50mg DTG is non-inferior to Truvada. Moreover, the predicted 90% preventive concentrations can guide release kinetics of currently developed DTG nano-formulations.

Diverse genetic error modes constrain large-scale bio-based production

A transition toward sustainable bio-based chemical production is important for green growth. However, productivity and yield frequently decrease as large-scale microbial fermentation progresses, commonly ascribed to phenotypic variation. Yet, given the high metabolic burden and toxicities, evolutionary processes may also constrain bio-based production. We experimentally simulate large-scale fermentation with mevalonic acid-producing Escherichia coli. By tracking growth rate and production, we uncover how populations fully sacrifice production to gain fitness within 70 generations. Using ultra-deep (>1000×) time-lapse sequencing of the pathway populations, we identify multiple recurring intra-pathway genetic error modes. This genetic heterogeneity is only detected using deep-sequencing and new population-level bioinformatics, suggesting that the problem is underestimated. A quantitative model explains the population dynamics based on enrichment of spontaneous mutant cells. We validate our model by tuning production load and escape rate of the production host and apply multiple orthogonal strategies for postponing genetically driven production declines.

The declining performance of scale-up bioreactor cultures is commonly attributed to phenotypic and physical heterogeneities. Here, the authors reveal multiple recurring intra-pathway error modes that limit engineered E. coli mevalonic acid production over time- and industrial-scale fermentations.

A stochastic and dynamical view of pluripotency in mouse embryonic stem cells

Pluripotent embryonic stem cells are of paramount importance for biomedical sciences because of their innate ability for self-renewal and differentiation into all major cell lines. The fateful decision to exit or remain in the pluripotent state is regulated by complex genetic regulatory networks. The rapid growth of single-cell sequencing data has greatly stimulated applications of statistical and machine learning methods for inferring topologies of pluripotency regulating genetic networks. The inferred network topologies, however, often only encode Boolean information while remaining silent about the roles of dynamics and molecular stochasticity inherent in gene expression. Herein we develop a framework for systematically extending Boolean-level network topologies into higher resolution models of networks which explicitly account for the promoter architectures and gene state switching dynamics. We show the framework to be useful for disentangling the various contributions that gene switching, external signaling, and network topology make to the global heterogeneity and dynamics of transcription factor populations. We find the pluripotent state of the network to be a steady state which is robust to global variations of gene switching rates which we argue are a good proxy for epigenetic states of individual promoters. The temporal dynamics of exiting the pluripotent state, on the other hand, is significantly influenced by the rates of genetic switching which makes cells more responsive to changes in extracellular signals.

Therapeutic target discovery using Boolean network attractors: improvements of kali

In a previous article, an algorithm for identifying therapeutic targets in Boolean networks modelling pathological mechanisms was introduced. In the present article, the improvements made on this algorithm, named kali, are described. These improvements are (i) the possibility to work on asynchronous Boolean networks, (ii) a finer assessment of therapeutic targets and (iii) the possibility to use multivalued logic. kali assumes that the attractors of a dynamical system, such as a Boolean network, are associated with the phenotypes of the modelled biological system. Given a logic-based model of pathological mechanisms, kali searches for therapeutic targets able to reduce the reachability of the attractors associated with pathological phenotypes, thus reducing their likeliness. kali is illustrated on an example network and used on a biological case study. The case study is a published logic-based model of bladder tumorigenesis from which kali returns consistent results. However, like any computational tool, kali can predict but cannot replace human expertise: it is a supporting tool for coping with the complexity of biological systems in the field of drug discovery.

A Multi-stage Representation of Cell Proliferation as a Markov Process

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.

LASSIE: simulating large-scale models of biochemical systems on GPUs

BACKGROUND

Mathematical modeling and in silico analysis are widely acknowledged as complementary tools to biological laboratory methods, to achieve a thorough understanding of emergent behaviors of cellular processes in both physiological and perturbed conditions. Though, the simulation of large-scale models-consisting in hundreds or thousands of reactions and molecular species-can rapidly overtake the capabilities of Central Processing Units (CPUs). The purpose of this work is to exploit alternative high-performance computing solutions, such as Graphics Processing Units (GPUs), to allow the investigation of these models at reduced computational costs.

RESULTS

LASSIE is a "black-box" GPU-accelerated deterministic simulator, specifically designed for large-scale models and not requiring any expertise in mathematical modeling, simulation algorithms or GPU programming. Given a reaction-based model of a cellular process, LASSIE automatically generates the corresponding system of Ordinary Differential Equations (ODEs), assuming mass-action kinetics. The numerical solution of the ODEs is obtained by automatically switching between the Runge-Kutta-Fehlberg method in the absence of stiffness, and the Backward Differentiation Formulae of first order in presence of stiffness. The computational performance of LASSIE are assessed using a set of randomly generated synthetic reaction-based models of increasing size, ranging from 64 to 8192 reactions and species, and compared to a CPU-implementation of the LSODA numerical integration algorithm.

CONCLUSIONS

LASSIE adopts a novel fine-grained parallelization strategy to distribute on the GPU cores all the calculations required to solve the system of ODEs. By virtue of this implementation, LASSIE achieves up to 92× speed-up with respect to LSODA, therefore reducing the running time from approximately 1 month down to 8 h to simulate models consisting in, for instance, four thousands of reactions and species. Notably, thanks to its smaller memory footprint, LASSIE is able to perform fast simulations of even larger models, whereby the tested CPU-implementation of LSODA failed to reach termination. LASSIE is therefore expected to make an important breakthrough in Systems Biology applications, for the execution of faster and in-depth computational analyses of large-scale models of complex biological systems.

Hybrid Markov-mass action law model for cell activation by rare binding events: Application to calcium induced vesicular release at neuronal synapses

Binding of molecules, ions or proteins to small target sites is a generic step of cell activation. This process relies on rare stochastic events where a particle located in a large bulk has to find small and often hidden targets. We present here a hybrid discrete-continuum model that takes into account a stochastic regime governed by rare events and a continuous regime in the bulk. The rare discrete binding events are modeled by a Markov chain for the encounter of small targets by few Brownian particles, for which the arrival time is Poissonian. The large ensemble of particles is described by mass action laws. We use this novel model to predict the time distribution of vesicular release at neuronal synapses. Vesicular release is triggered by the binding of few calcium ions that can originate either from the synaptic bulk or from the entry through calcium channels. We report here that the distribution of release time is bimodal although it is triggered by a single fast action potential. While the first peak follows a stimulation, the second corresponds to the random arrival over much longer time of ions located in the synaptic terminal to small binding vesicular targets. To conclude, the present multiscale stochastic modeling approach allows studying cellular events based on integrating discrete molecular events over several time scales.

Life sciences today and tomorrow: emerging biotechnologies

The purpose of this review is to survey current, emerging and predicted future biotechnologies which are impacting, or are likely to impact in the future on the life sciences, with a projection for the coming 20 years. This review is intended to discuss current and future technical strategies, and to explore areas of potential growth during the foreseeable future. Information technology approaches have been employed to gather and collate data. Twelve broad categories of biotechnology have been identified which are currently impacting the life sciences and will continue to do so. In some cases, technology areas are being pushed forward by the requirement to deal with contemporary questions such as the need to address the emergence of anti-microbial resistance. In other cases, the biotechnology application is made feasible by advances in allied fields in biophysics (e.g. biosensing) and biochemistry (e.g. bio-imaging). In all cases, the biotechnologies are underpinned by the rapidly advancing fields of information systems, electronic communications and the World Wide Web together with developments in computing power and the capacity to handle extensive biological data. A rationale and narrative is given for the identification of each technology as a growth area. These technologies have been categorized by major applications, and are discussed further. This review highlights: Biotechnology has far-reaching applications which impinge on every aspect of human existence. The applications of biotechnology are currently wide ranging and will become even more diverse in the future. Access to supercomputing facilities and the ability to manipulate large, complex biological datasets, will significantly enhance knowledge and biotechnological development.

Novel Approaches Reveal that Toxoplasma gondii Bradyzoites within Tissue Cysts Are Dynamic and Replicating Entities In Vivo

Despite their critical role in chronic toxoplasmosis, the biology of Toxoplasma gondii bradyzoites is poorly understood. In an attempt to address this gap, we optimized approaches to purify tissue cysts and analyzed the replicative potential of bradyzoites within these cysts. In order to quantify individual bradyzoites within tissue cysts, we have developed imaging software, BradyCount 1.0, that allows the rapid establishment of bradyzoite burdens within imaged optical sections of purified tissue cysts. While in general larger tissue cysts contain more bradyzoites, their relative “occupancy” was typically lower than that of smaller cysts, resulting in a lower packing density. The packing density permits a direct measure of how bradyzoites develop within cysts, allowing for comparisons across progression of the chronic phase. In order to capture bradyzoite endodyogeny, we exploited the differential intensity of TgIMC3, an inner membrane complex protein that intensely labels newly formed/forming daughters within bradyzoites and decays over time in the absence of further division. To our surprise, we were able to capture not only sporadic and asynchronous division but also synchronous replication of all bradyzoites within mature tissue cysts. Furthermore, the time-dependent decay of TgIMC3 intensity was exploited to gain insights into the temporal patterns of bradyzoite replication in vivo. Despite the fact that bradyzoites are considered replicatively dormant, we find evidence for cyclical, episodic bradyzoite growth within tissue cysts in vivo. These findings directly challenge the prevailing notion of bradyzoites as dormant nonreplicative entities in chronic toxoplasmosis and have implications on our understanding of this enigmatic and clinically important life cycle stage.

The protozoan Toxoplasma gondii establishes a lifelong chronic infection mediated by the bradyzoite form of the parasite within tissue cysts. Technical challenges have limited even the most basic studies on bradyzoites and the tissue cysts in vivo. Bradyzoites, which are viewed as dormant, poorly replicating or nonreplicating entities, were found to be surprisingly active, exhibiting not only the capacity for growth but also previously unrecognized patterns of replication that point to their being considerably more dynamic than previously imagined. These newly revealed properties force us to reexamine the most basic questions regarding bradyzoite biology and the progression of the chronic phase of toxoplasmosis. By developing new tools and approaches to study the chronic phase at the level of bradyzoites, we expose new avenues to tackle both drug development and a better understanding of events that may lead to reactivated symptomatic disease.

A Computational Framework for Bioimaging Simulation

Using bioimaging technology, biologists have attempted to identify and document analytical interpretations that underlie biological phenomena in biological cells. Theoretical biology aims at distilling those interpretations into knowledge in the mathematical form of biochemical reaction networks and understanding how higher level functions emerge from the combined action of biomolecules. However, there still remain formidable challenges in bridging the gap between bioimaging and mathematical modeling. Generally, measurements using fluorescence microscopy systems are influenced by systematic effects that arise from stochastic nature of biological cells, the imaging apparatus, and optical physics. Such systematic effects are always present in all bioimaging systems and hinder quantitative comparison between the cell model and bioimages. Computational tools for such a comparison are still unavailable. Thus, in this work, we present a computational framework for handling the parameters of the cell models and the optical physics governing bioimaging systems. Simulation using this framework can generate digital images of cell simulation results after accounting for the systematic effects. We then demonstrate that such a framework enables comparison at the level of photon-counting units.

Fractal model of anomalous diffusion

An equation of motion is derived from fractal analysis of the Brownian particle trajectory in which the asymptotic fractal dimension of the trajectory has a required value. The formula makes it possible to calculate the time dependence of the mean square displacement for both short and long periods when the molecule diffuses anomalously. The anomalous diffusion which occurs after long periods is characterized by two variables, the transport coefficient and the anomalous diffusion exponent. An explicit formula is derived for the transport coefficient, which is related to the diffusion constant, as dependent on the Brownian step time, and the anomalous diffusion exponent. The model makes it possible to deduce anomalous diffusion properties from experimental data obtained even for short time periods and to estimate the transport coefficient in systems for which the diffusion behavior has been investigated. The results were confirmed for both sub and super-diffusion.