Search

We present StochSS: Stochastic Simulation as a Service, an integrated development environment for modeling and simulation of both deterministic and discrete stochastic biochemical systems in up to three dimensions. An easy to use graphical user interface enables researchers to quickly develop and simulate a biological model on a desktop or laptop, which can then be expanded to incorporate increasing levels of complexity. StochSS features state-of-the-art simulation engines. As the demand for computational power increases, StochSS can seamlessly scale computing resources in the cloud. In addition, StochSS can be deployed as a multi-user software environment where collaborators share computational resources and exchange models via a public model repository. We demonstrate the capabilities and ease of use of StochSS with an example of model development and simulation at increasing levels of complexity.

Cells polarize their movement or growth toward external directional cues in many different contexts. For example, budding yeast cells grow toward potential mating partners in response to pheromone gradients. Directed growth is controlled by polarity factors that assemble into clusters at the cell membrane. The clusters assemble, disassemble, and move between different regions of the membrane before eventually forming a stable polarity site directed toward the pheromone source. Pathways that regulate clustering have been identified but the molecular mechanisms that regulate cluster mobility are not well understood. To gain insight into the contribution of chemical noise to cluster behavior we simulated clustering using the reaction-diffusion master equation (RDME) framework to account for molecular-level fluctuations. RDME simulations are a computationally efficient approximation, but their results can diverge from the underlying microscopic dynamics. We implemented novel concentration-dependent rate constants that improved the accuracy of RDME-based simulations, allowing us to efficiently investigate how cluster dynamics might be regulated. Molecular noise was effective in relocating clusters when the clusters contained low numbers of limiting polarity factors, and when Cdc42, the central polarity regulator, exhibited short dwell times at the polarity site. Cluster stabilization occurred when abundances or binding rates were altered to either lengthen dwell times or increase the number of polarity molecules in the cluster. We validated key results using full 3D particle-based simulations. Understanding the mechanisms cells use to regulate the dynamics of polarity clusters should provide insights into how cells dynamically track external directional cues.

During the last decade, a consensus has emerged that the stochastic triggering of an excitable system drives pseudopod formation and subsequent migration of amoeboid cells. The presence of chemoattractant stimuli alters the threshold for triggering this activity and can bias the direction of migration. Though noise plays an important role in these behaviors, mathematical models have typically ignored its origin and merely introduced it as an external signal into a series of reaction-diffusion equations. Here we consider a more realistic description based on a reaction-diffusion master equation formalism to implement these networks. In this scheme, noise arises naturally from a stochastic description of the various reaction and diffusion terms. Working on a three-dimensional geometry in which separate compartments are divided into a tetrahedral mesh, we implement a modular description of the system, consisting of G-protein coupled receptor signaling (GPCR), a local excitation-global inhibition mechanism (LEGI), and signal transduction excitable network (STEN). Our models implement detailed biochemical descriptions whenever this information is available, such as in the GPCR and G-protein interactions. In contrast, where the biochemical entities are less certain, such as the LEGI mechanism, we consider various possible schemes and highlight the differences between them. Our simulations show that even when the LEGI mechanism displays perfect adaptation in terms of the mean level of proteins, the variance shows a dose-dependence. This differs between the various models considered, suggesting a possible means for determining experimentally among the various potential networks. Overall, our simulations recreate temporal and spatial patterns observed experimentally in both wild-type and perturbed cells, providing further evidence for the excitable system paradigm. Moreover, because of the overall importance and ubiquity of the modules we consider, including GPCR signaling and adaptation, our results will be of interest beyond the field of directed migration.

Skeletal muscle contraction is initiated when an action potential triggers the release of Ca2+ into the sarcomere in a process referred to as excitation-contraction coupling. The speed and scale of this process makes direct observation very challenging and invasive. To determine how the concentration of Ca2+ changes within the myofibril during a single activation, several simulation models have been developed. These models follow a common pattern; divide the half sarcomere into a series of compartments, then use ordinary differential equations to solve reactions occurring within and between the compartments. To further develop this type of simulation, we have created a realistic structural model of a skeletal muscle myofibrillar half-sarcomere using MCell software that incorporates the myofilament lattice structure. Using this simulation model, we were successful in reproducing the averaged calcium transient during a single activation consistent with both the experimental and previous simulation results. In addition, our simulation demonstrated that the inclusion of the myofilament lattice within our model produced an asymmetric distribution of Ca2+, with more Ca2+ accumulating near the Z-disk and less Ca2+ reaching the m-line. This asymmetric distribution of Ca2+ is also apparent when we examine how the Ca2+ are bound to the troponin-C proteins along the actin filaments. Our simulation model also allowed us to produce advanced visualizations of this process, including two simulation animations, allowing us to view Ca2+ release, diffusion, binding and uptake within the myofibrillar half-sarcomere.

Calcium/calmodulin-dependent protein kinase II (CaMKII) holoenzymes play a critical role in decoding Ca2+ signals in neurons. Understanding how this occurs has been the focus of numerous studies including many that use models. However, CaMKII is notoriously difficult to simulate in detail because of its multi-subunit nature, which causes a combinatorial explosion in the number of species that must be modeled. To study the Ca2+-calmodulin-CaMKII reaction network with detailed kinetics while including the effect of diffusion, we have customized an existing stochastic particle-based simulator, Smoldyn, to manage the problem of combinatorial explosion. With this new method, spatial and temporal aspects of the signaling network can be studied without compromising biochemical details. We used this new method to examine how calmodulin molecules, both partially loaded and fully loaded with Ca2+, choose pathways to interact with and activate CaMKII under various Ca2+ input conditions. We found that the dependence of CaMKII phosphorylation on Ca2+ signal frequency is intrinsic to the network kinetics and the activation pattern can be modulated by the relative amount of Ca2+ to calmodulin and by the rate of Ca2+ diffusion. Depending on whether Ca2+ influx is saturating or not, calmodulin molecules could choose different routes within the network to activate CaMKII subunits, resulting in different frequency dependence patterns. In addition, the size of the holoenzyme produces a subtle effect on CaMKII activation. The more extended the subunits are organized, the easier for calmodulin molecules to access and activate the subunits. The findings suggest that particular intracellular environmental factors such as crowding and calmodulin availability can play an important role in decoding Ca2+ signals and can give rise to distinct CaMKII activation patterns in dendritic spines, Ca2+ channel nanodomains and cytoplasm.

We present StochSS: Stochastic Simulation as a Service, an integrated development environment for modeling and simulation of both deterministic and discrete stochastic biochemical systems in up to three dimensions. An easy to use graphical user interface enables researchers to quickly develop and simulate a biological model on a desktop or laptop, which can then be expanded to incorporate increasing levels of complexity. StochSS features state-of-the-art simulation engines. As the demand for computational power increases, StochSS can seamlessly scale computing resources in the cloud. In addition, StochSS can be deployed as a multi-user software environment where collaborators share computational resources and exchange models via a public model repository. We demonstrate the capabilities and ease of use of StochSS with an example of model development and simulation at increasing levels of complexity.

Cells polarize their movement or growth toward external directional cues in many different contexts. For example, budding yeast cells grow toward potential mating partners in response to pheromone gradients. Directed growth is controlled by polarity factors that assemble into clusters at the cell membrane. The clusters assemble, disassemble, and move between different regions of the membrane before eventually forming a stable polarity site directed toward the pheromone source. Pathways that regulate clustering have been identified but the molecular mechanisms that regulate cluster mobility are not well understood. To gain insight into the contribution of chemical noise to cluster behavior we simulated clustering using the reaction-diffusion master equation (RDME) framework to account for molecular-level fluctuations. RDME simulations are a computationally efficient approximation, but their results can diverge from the underlying microscopic dynamics. We implemented novel concentration-dependent rate constants that improved the accuracy of RDME-based simulations, allowing us to efficiently investigate how cluster dynamics might be regulated. Molecular noise was effective in relocating clusters when the clusters contained low numbers of limiting polarity factors, and when Cdc42, the central polarity regulator, exhibited short dwell times at the polarity site. Cluster stabilization occurred when abundances or binding rates were altered to either lengthen dwell times or increase the number of polarity molecules in the cluster. We validated key results using full 3D particle-based simulations. Understanding the mechanisms cells use to regulate the dynamics of polarity clusters should provide insights into how cells dynamically track external directional cues., Competing Interests: The authors have declared that no competing interests exist.

During the last decade, a consensus has emerged that the stochastic triggering of an excitable system drives pseudopod formation and subsequent migration of amoeboid cells. The presence of chemoattractant stimuli alters the threshold for triggering this activity and can bias the direction of migration. Though noise plays an important role in these behaviors, mathematical models have typically ignored its origin and merely introduced it as an external signal into a series of reaction-diffusion equations. Here we consider a more realistic description based on a reaction-diffusion master equation formalism to implement these networks. In this scheme, noise arises naturally from a stochastic description of the various reaction and diffusion terms. Working on a three-dimensional geometry in which separate compartments are divided into a tetrahedral mesh, we implement a modular description of the system, consisting of G-protein coupled receptor signaling (GPCR), a local excitation-global inhibition mechanism (LEGI), and signal transduction excitable network (STEN). Our models implement detailed biochemical descriptions whenever this information is available, such as in the GPCR and G-protein interactions. In contrast, where the biochemical entities are less certain, such as the LEGI mechanism, we consider various possible schemes and highlight the differences between them. Our simulations show that even when the LEGI mechanism displays perfect adaptation in terms of the mean level of proteins, the variance shows a dose-dependence. This differs between the various models considered, suggesting a possible means for determining experimentally among the various potential networks. Overall, our simulations recreate temporal and spatial patterns observed experimentally in both wild-type and perturbed cells, providing further evidence for the excitable system paradigm. Moreover, because of the overall importance and ubiquity of the modules we consider, including GPCR signaling and adaptation, our results will be of interest beyond the field of directed migration., Competing Interests: The authors have declared that no competing interests exist.

Simulation of biomolecular networks is now indispensable for studying biological systems, from small reaction networks to large ensembles of cells. Here we present a novel approach for stochastic simulation of networks embedded in the dynamic environment of the cell and its surroundings. We thus sample trajectories of the stochastic process described by the chemical master equation with time-varying propensities. A comparative analysis shows that existing approaches can either fail dramatically, or else can impose impractical computational burdens due to numerical integration of reaction propensities, especially when cell ensembles are studied. Here we introduce the Extrande method which, given a simulated time course of dynamic network inputs, provides a conditionally exact and several orders-of-magnitude faster simulation solution. The new approach makes it feasible to demonstrate-using decision-making by a large population of quorum sensing bacteria-that robustness to fluctuations from upstream signaling places strong constraints on the design of networks determining cell fate. Our approach has the potential to significantly advance both understanding of molecular systems biology and design of synthetic circuits.

Skeletal muscle contraction is initiated when an action potential triggers the release of Ca2+ into the sarcomere in a process referred to as excitation-contraction coupling. The speed and scale of this process makes direct observation very challenging and invasive. To determine how the concentration of Ca2+ changes within the myofibril during a single activation, several simulation models have been developed. These models follow a common pattern; divide the half sarcomere into a series of compartments, then use ordinary differential equations to solve reactions occurring within and between the compartments. To further develop this type of simulation, we have created a realistic structural model of a skeletal muscle myofibrillar half-sarcomere using MCell software that incorporates the myofilament lattice structure. Using this simulation model, we were successful in reproducing the averaged calcium transient during a single activation consistent with both the experimental and previous simulation results. In addition, our simulation demonstrated that the inclusion of the myofilament lattice within our model produced an asymmetric distribution of Ca2+, with more Ca2+ accumulating near the Z-disk and less Ca2+ reaching the m-line. This asymmetric distribution of Ca2+ is also apparent when we examine how the Ca2+ are bound to the troponin-C proteins along the actin filaments. Our simulation model also allowed us to produce advanced visualizations of this process, including two simulation animations, allowing us to view Ca2+ release, diffusion, binding and uptake within the myofibrillar half-sarcomere., Competing Interests: The authors have declared that no competing interests exist.

Calcium/calmodulin-dependent protein kinase II (CaMKII) holoenzymes play a critical role in decoding Ca2+ signals in neurons. Understanding how this occurs has been the focus of numerous studies including many that use models. However, CaMKII is notoriously difficult to simulate in detail because of its multi-subunit nature, which causes a combinatorial explosion in the number of species that must be modeled. To study the Ca2+-calmodulin-CaMKII reaction network with detailed kinetics while including the effect of diffusion, we have customized an existing stochastic particle-based simulator, Smoldyn, to manage the problem of combinatorial explosion. With this new method, spatial and temporal aspects of the signaling network can be studied without compromising biochemical details. We used this new method to examine how calmodulin molecules, both partially loaded and fully loaded with Ca2+, choose pathways to interact with and activate CaMKII under various Ca2+ input conditions. We found that the dependence of CaMKII phosphorylation on Ca2+ signal frequency is intrinsic to the network kinetics and the activation pattern can be modulated by the relative amount of Ca2+ to calmodulin and by the rate of Ca2+ diffusion. Depending on whether Ca2+ influx is saturating or not, calmodulin molecules could choose different routes within the network to activate CaMKII subunits, resulting in different frequency dependence patterns. In addition, the size of the holoenzyme produces a subtle effect on CaMKII activation. The more extended the subunits are organized, the easier for calmodulin molecules to access and activate the subunits. The findings suggest that particular intracellular environmental factors such as crowding and calmodulin availability can play an important role in decoding Ca2+ signals and can give rise to distinct CaMKII activation patterns in dendritic spines, Ca2+ channel nanodomains and cytoplasm.

Simulation of biomolecular networks is now indispensable for studying biological systems, from small reaction networks to large ensembles of cells. Here we present a novel approach for stochastic simulation of networks embedded in the dynamic environment of the cell and its surroundings. We thus sample trajectories of the stochastic process described by the chemical master equation with time-varying propensities. A comparative analysis shows that existing approaches can either fail dramatically, or else can impose impractical computational burdens due to numerical integration of reaction propensities, especially when cell ensembles are studied. Here we introduce the Extrande method which, given a simulated time course of dynamic network inputs, provides a conditionally exact and several orders-of-magnitude faster simulation solution. The new approach makes it feasible to demonstrate—using decision-making by a large population of quorum sensing bacteria—that robustness to fluctuations from upstream signaling places strong constraints on the design of networks determining cell fate. Our approach has the potential to significantly advance both understanding of molecular systems biology and design of synthetic circuits. [ABSTRACT FROM AUTHOR]

Although the length and constituting sequences for pericentromeric repeats are highly variable across eukaryotes, the presence of multiple pericentromeric repeats is one of the conserved features of the eukaryotic chromosomes. Pericentromeric heterochromatin is often misregulated in human diseases, with the expansion of pericentromeric repeats in human solid cancers. In this article, we have developed a mathematical model of the RNAi-dependent methylation of H3K9 in the pericentromeric region of fission yeast. Our model, which takes copy number as an explicit parameter, predicts that the pericentromere is silenced only if there are many copies of repeats. It becomes bistable or desilenced if the copy number of repeats is reduced. This suggests that the copy number of pericentromeric repeats alone can determine the fate of heterochromatin silencing in fission yeast. Through sensitivity analysis, we identified parameters that favor bistability and desilencing. Stochastic simulation shows that faster cell division and noise favor the desilenced state. These results show the unexpected role of pericentromeric repeat copy number in gene silencing and provide a quantitative basis for how the copy number allows or protects repetitive and unique parts of the genome from heterochromatin silencing, respectively. Author summary: Pericentromeric repeats vary in length and sequences, but their presence is a conserved feature of eukaryotes. This suggests that the repetitive nature of pericentromeric sequences is a functional feature of centromeres, which is under selective pressure. Here we developed a quantitative model for gene silencing at the fission yeast pericentromeric repeats. Our model is one of the first models which incorporates the copy number of pericentromeric repeats and predicts that the number of repeats can solely govern the dynamics of pericentromeric gene silencing. Our results suggest that the repeat copy number is critical for gene silencing, and copy-number-dependent silencing is an effective strategy to repress the repetitive part of the genome and protect the unique part of the genome from heterochromatin silencing. [ABSTRACT FROM AUTHOR]

Pre-exposure prophylaxis (PrEP) is an important pillar to prevent HIV transmission. Because of experimental and clinical shortcomings, mathematical models that integrate pharmacological, viral- and host factors are frequently used to quantify clinical efficacy of PrEP. Stochastic simulations of these models provides sample statistics from which the clinical efficacy is approximated. However, many stochastic simulations are needed to reduce the associated sampling error. To remedy the shortcomings of stochastic simulation, we developed a numerical method that allows predicting the efficacy of arbitrary prophylactic regimen directly from a viral dynamics model, without sampling. We apply the method to various hypothetical dolutegravir (DTG) prophylaxis scenarios. The approach is verified against state-of-the-art stochastic simulation. While the method is more accurate than stochastic simulation, it is superior in terms of computational performance. For example, a continuous 6-month prophylactic profile is computed within a few seconds on a laptop computer. The method's computational performance, therefore, substantially expands the horizon of feasible analysis in the context of PrEP, and possibly other applications. Author summary: Pre-exposure prophylaxis (PrEP) is an important tool to prevent HIV transmission. However, experimental identification of parameters that determine prophylactic efficacy is extremely difficult. Clues about these parameters could prove essential for the design of next-generation PrEP compounds. Integrative mathematical models can fill this void: Based on stochastic simulation, a sample statistic can be generated, from which the prophylactic efficacy is estimated. However, for this sample statistic to be accurate, many simulations need to be performed. Here, we introduce a numerical method to directly compute the prophylactic efficacy from a viral dynamics model, without the need for sampling. Based on several examples with dolutegravir (DTG) -based short- and long-term PrEP, as well as post-exposure prophylaxis we demonstrate the correctness of the new method and its outstanding computational performance. Due to the method's computational performance, a number of analyses, including formal sensitivity analysis, are becoming feasible with the proposed method. [ABSTRACT FROM AUTHOR]

Patients with coronavirus disease 2019 (COVID-19) often exhibit diverse disease progressions associated with various infectious ability, symptoms, and clinical treatments. To systematically and thoroughly understand the heterogeneous progression of COVID-19, we developed a multi-scale computational model to quantitatively understand the heterogeneous progression of COVID-19 patients infected with severe acute respiratory syndrome (SARS)-like coronavirus (SARS-CoV-2). The model consists of intracellular viral dynamics, multicellular infection process, and immune responses, and was formulated using a combination of differential equations and stochastic modeling. By integrating multi-source clinical data with model analysis, we quantified individual heterogeneity using two indexes, i.e., the ratio of infected cells and incubation period. Specifically, our simulations revealed that increasing the host antiviral state or virus induced type I interferon (IFN) production rate can prolong the incubation period and postpone the transition from asymptomatic to symptomatic outcomes. We further identified the threshold dynamics of T cell exhaustion in the transition between mild-moderate and severe symptoms, and that patients with severe symptoms exhibited a lack of naïve T cells at a late stage. In addition, we quantified the efficacy of treating COVID-19 patients and investigated the effects of various therapeutic strategies. Simulations results suggested that single antiviral therapy is sufficient for moderate patients, while combination therapies and prevention of T cell exhaustion are needed for severe patients. These results highlight the critical roles of IFN and T cell responses in regulating the stage transition during COVID-19 progression. Our study reveals a quantitative relationship underpinning the heterogeneity of transition stage during COVID-19 progression and can provide a potential guidance for personalized therapy in COVID-19 patients. Author summary: Coronavirus disease 2019 (COVID-19) is currently destroying both lives and economies. However, patients infected with severe acute respiratory syndrome (SARS)-like coronavirus (SARS-CoV-2) usually present heterogeneous and complicated progressions, such as different incubation periods (short and long), symptoms (asymptomatic and symptomatic) and severity (mild-moderate and severe). Currently, various clinical data and experimental data are available from different countries, which has great significance for integrating different types of data to comprehensively understand the diverse disease progression in COVID-19 patients and guide individual treatment strategies. Here, we developed a multi-scale computational model to describe the dynamical process of patients infected with SARS-CoV-2, including intracellular viral dynamics, multicellular infection process, and immune responses. By combining data integration, stochastic simulation and quantitative analysis based on the multi-scale mathematical model, we addressed an important question regarding how IFN response and T cell exhaustion quantitatively affect heterogeneous progression in patients with respect to incubation periods, symptoms and severity. Furthermore, the efficacy of various therapeutic strategies for treating COVID-19 patients with different severity degrees was evaluated and validated. The computational framework in this study can also be extended to explore the dynamical process of other coronavirus infections. [ABSTRACT FROM AUTHOR]

Recombinant protein production is a key process in generating proteins of interest in the pharmaceutical industry and biomedical research. However, about 50% of recombinant proteins fail to be expressed in a variety of host cells. Here we show that the accessibility of translation initiation sites modelled using the mRNA base-unpairing across the Boltzmann's ensemble significantly outperforms alternative features. This approach accurately predicts the successes or failures of expression experiments, which utilised Escherichia coli cells to express 11,430 recombinant proteins from over 189 diverse species. On this basis, we develop TIsigner that uses simulated annealing to modify up to the first nine codons of mRNAs with synonymous substitutions. We show that accessibility captures the key propensity beyond the target region (initiation sites in this case), as a modest number of synonymous changes is sufficient to tune the recombinant protein expression levels. We build a stochastic simulation model and show that higher accessibility leads to higher protein production and slower cell growth, supporting the idea of protein cost, where cell growth is constrained by protein circuits during over expression. Author summary: Recombinant proteins are widely used as therapeutics, such as vaccines, monoclonal antibodies, hormones and enzymes. However, the success rate of recombinant protein production is about 50%. To address this problem, we propose optimising the unpairing propensities of nucleotides around translation initiation sites using a thermodynamic quantity called mRNA accessibility. Our study shows that this method is generalisable across prokaryotic and eukaryotic expression hosts. Importantly, we validated this method using laboratory experiments and computational modelling. Furthermore, we propose a low cost technique to tune protein expression by engineering minimal changes to genes of interest through our web application (https://tisigner.com/tisigner). [ABSTRACT FROM AUTHOR]

Stochastic simulation has been a powerful tool for studying the dynamics of gene regulatory networks, particularly in terms of understanding how cell-phenotype stability and fate-transitions are impacted by noisy gene expression. However, gene networks often have dynamics characterized by multiple attractors. Stochastic simulation is often inefficient for such systems, because most of the simulation time is spent waiting for rare, barrier-crossing events to occur. We present a rare-event simulation-based method for computing epigenetic landscapes and phenotype-transitions in metastable gene networks. Our computational pipeline was inspired by studies of metastability and barrier-crossing in protein folding, and provides an automated means of computing and visualizing essential stationary and dynamic information that is generally inaccessible to conventional simulation. Applied to a network model of pluripotency in Embryonic Stem Cells, our simulations revealed rare phenotypes and approximately Markovian transitions among phenotype-states, occurring with a broad range of timescales. The relative probabilities of phenotypes and the transition paths linking pluripotency and differentiation are sensitive to global kinetic parameters governing transcription factor-DNA binding kinetics. Our approach significantly expands the capability of stochastic simulation to investigate gene regulatory network dynamics, which may help guide rational cell reprogramming strategies. Our approach is also generalizable to other types of molecular networks and stochastic dynamics frameworks. [ABSTRACT FROM AUTHOR]

There is an abundance of malaria genetic data being collected from the field, yet using these data to understand the drivers of regional epidemiology remains a challenge. A key issue is the lack of models that relate parasite genetic diversity to epidemiological parameters. Classical models in population genetics characterize changes in genetic diversity in relation to demographic parameters, but fail to account for the unique features of the malaria life cycle. In contrast, epidemiological models, such as the Ross-Macdonald model, capture malaria transmission dynamics but do not consider genetics. Here, we have developed an integrated model encompassing both parasite evolution and regional epidemiology. We achieve this by combining the Ross-Macdonald model with an intra-host continuous-time Moran model, thus explicitly representing the evolution of individual parasite genomes in a traditional epidemiological framework. Implemented as a stochastic simulation, we use the model to explore relationships between measures of parasite genetic diversity and parasite prevalence, a widely-used metric of transmission intensity. First, we explore how varying parasite prevalence influences genetic diversity at equilibrium. We find that multiple genetic diversity statistics are correlated with prevalence, but the strength of the relationships depends on whether variation in prevalence is driven by host- or vector-related factors. Next, we assess the responsiveness of a variety of statistics to malaria control interventions, finding that those related to mixed infections respond quickly (∼months) whereas other statistics, such as nucleotide diversity, may take decades to respond. These findings provide insights into the opportunities and challenges associated with using genetic data to monitor malaria epidemiology. Author Summary: Knowledge of how the prevalence of P.falciparum malaria varies, either between regions or through time, is critical to the operation of malaria control programs. Yet obtaining this information through traditional methods presents many challenges. Parasite genetic data is increasingly accessible, and may provide an alternative means to estimate P.falciparum prevalence in the field. However, our understanding of how the genetic diversity of parasite populations relates to prevalence is limited, and suitable models to guide our understanding are largely lacking. Here, we merge two classical models—the Ross-Macondald and the Moran—to produce a framework in which the relationships between parasite genetic diversity and prevalence can be explored. We find that several genetic diversity statistics are correlated with prevalence, although to differing degrees, and over different time scales. Overall, statistics related to mixed infection are robustly and rapidly responsive to changes in prevalence, suggesting they may be a useful focal point for the development of malaria surveillance methods that harness genetic data. [ABSTRACT FROM AUTHOR]

Chemical reaction networks are ubiquitous in biology, and their dynamics is fundamentally stochastic. Here, we present the software library pSSAlib, which provides a complete and concise implementation of the most efficient partial-propensity methods for simulating exact stochastic chemical kinetics. pSSAlib can import models encoded in Systems Biology Markup Language, supports time delays in chemical reactions, and stochastic spatiotemporal reaction-diffusion systems. It also provides tools for statistical analysis of simulation results and supports multiple output formats. It has previously been used for studies of biochemical reaction pathways and to benchmark other stochastic simulation methods. Here, we describe pSSAlib in detail and apply it to a new model of the endocytic pathway in eukaryotic cells, leading to the discovery of a stochastic counterpart of the cut-out switch motif underlying early-to-late endosome conversion. pSSAlib is provided as a stand-alone command-line tool and as a developer API. We also provide a plug-in for the SBMLToolbox. The open-source code and pre-packaged installers are freely available from . [ABSTRACT FROM AUTHOR]

Accurately estimating the effective reproduction number (Rt) of a circulating pathogen is a fundamental challenge in the study of infectious disease. The fields of epidemiology and pathogen phylodynamics both share this goal, but to date, methodologies and data employed by each remain largely distinct. Here we present EpiFusion: a joint approach that can be used to harness the complementary strengths of each field to improve estimation of outbreak dynamics for large and poorly sampled epidemics, such as arboviral or respiratory virus outbreaks, and validate it for retrospective analysis. We propose a model of Rt that estimates outbreak trajectories conditional upon both phylodynamic (time-scaled trees estimated from genetic sequences) and epidemiological (case incidence) data. We simulate stochastic outbreak trajectories that are weighted according to epidemiological and phylodynamic observation models and fit using particle Markov Chain Monte Carlo. To assess performance, we test EpiFusion on simulated outbreaks in which transmission and/or surveillance rapidly changes and find that using EpiFusion to combine epidemiological and phylodynamic data maintains accuracy and increases certainty in trajectory and Rt estimates, compared to when each data type is used alone. We benchmark EpiFusion's performance against existing methods to estimate Rt and demonstrate advances in speed and accuracy. Importantly, our approach scales efficiently with dataset size. Finally, we apply our model to estimate Rt during the 2014 Ebola outbreak in Sierra Leone. EpiFusion is designed to accommodate future extensions that will improve its utility, such as explicitly modelling population structure, accommodations for phylogenetic uncertainty, and the ability to weight the contributions of genomic or case incidence to the inference. Author summary: Understanding infectious disease spread is fundamental to protecting public health, but can be challenging as disease spread is a phenomenon that cannot be directly observed. So, epidemiologists use data in conjunction with mathematical models to estimate disease dynamics. Often, combinations of different models and data can be used to answer the same questions–for example 'traditional' epidemiology commonly uses case incidence data (the number of people who have tested positive for a disease during a certain time period) whereas phylodynamic models use pathogen genomic sequence data and our knowledge of the way their genomes evolve to model disease population dynamics. Each of these approaches have strengths and limitations, and data of each type can be sparse or biased, particularly during rapidly developing outbreaks or in countries with poor pathogen surveillance. An increasing number of approaches attempt to fix this problem by incorporating diverse concepts and data types together in their models. We aim to contribute to this movement by introducing EpiFusion, a modelling framework that improves the efficiency and precision at which we can monitor important changes in pathogen transmission (specifically, in the effective reproduction number). EpiFusion uses particle filtering to simulate epidemic trajectories over time and weight their likelihood according to both case incidence data and a phylogenetic tree using separate observation models, resulting in the inference of trajectories in agreement with both sets of data. Improvements in our ability to accurately and confidently model pathogen spread help us to respond to infectious disease outbreaks and improve public health. [ABSTRACT FROM AUTHOR]

The fate choices of stem cells between self-renewal and differentiation are often tightly regulated by juxtacrine (cell-cell contact) signalling. Here, we assess how the interplay between cell division, cell fate choices, and juxtacrine signalling can affect the macroscopic ordering of cell types in self-renewing epithelial sheets, by studying a simple spatial cell fate model with cells being arranged on a 2D lattice. We show in this model that if cells commit to their fate directly upon cell division, macroscopic patches of cells of the same type emerge, if at least a small proportion of divisions are symmetric, except if signalling interactions are laterally inhibiting. In contrast, if cells are first 'licensed' to differentiate, yet retaining the possibility to return to their naive state, macroscopic order only emerges if the signalling strength exceeds a critical threshold: if then the signalling interactions are laterally inducing, macroscopic patches emerge as well. Lateral inhibition, on the other hand, can in that case generate periodic patterns of alternating cell types (checkerboard pattern), yet only if the proportion of symmetric divisions is sufficiently low. These results can be understood theoretically by an analogy to phase transitions in spin systems known from statistical physics. Author summary: A fundamental question in stem cell biology is how a cell's choice to differentiate or not (cell fate choice), is regulated by communication with other cells in a tissue, and whether these choices are a one-way path or to some degree reversible. However, measuring this in living animals is very difficult and often impossible, since this requires to make videos of cells inside the body with a microscope. Here, we employ a simple mathematical model for the fate choices of stem cells when they are regulated by communication with nearby cells in the tissue. We show that different means of cell fate choice and cell communication can lead to qualitatively different macroscopic features of the spatial arrangement of cell types: large patches, checkerboard patterns, or randomly disordered distributions, depending on the character of cell communication, and whether cell fate is committed at cell division or reversible. Our analysis therefore shows that those aspects of stem cell activity, which are otherwise difficult to measure, can be distinguished by observing spatial arrangements of cell types. [ABSTRACT FROM AUTHOR]

Symbolic regression with polynomial neural networks and polynomial neural ordinary differential equations (ODEs) are two recent and powerful approaches for equation recovery of many science and engineering problems. However, these methods provide point estimates for the model parameters and are currently unable to accommodate noisy data. We address this challenge by developing and validating the following Bayesian inference methods: the Laplace approximation, Markov Chain Monte Carlo (MCMC) sampling methods, and variational inference. We have found the Laplace approximation to be the best method for this class of problems. Our work can be easily extended to the broader class of symbolic neural networks to which the polynomial neural network belongs. Author summary: Polynomial neural ordinary differential equations (ODEs) are a recent approach for symbolic regression of dynamical systems governed by polynomials. However, they are limited in that they provide maximum likelihood point estimates of the model parameters. The domain expert using system identification often desires a specified level of confidence or range of parameter values that best fit the data. In this work, we use Bayesian inference to provide posterior probability distributions of the parameters in polynomial neural ODEs. To date, there are no studies that attempt to identify the best Bayesian inference method for neural ODEs and symbolic neural ODEs. To address this need, we explore and compare three different approaches for estimating the posterior distributions of weights and biases of the polynomial neural network: the Laplace approximation, Markov Chain Monte Carlo (MCMC) sampling, and variational inference. We have found the Laplace approximation to be the best method for this class of problems. We have also developed lightweight JAX code to estimate posterior probability distributions using the Laplace approximation. [ABSTRACT FROM AUTHOR]

We discuss the invasion of the Omicron BA.1 variant into England as a paradigm for real-time model fitting and projection. Here we use a mixture of simple SIR-type models, analysis of the early data and a more complex age-structure model fit to the outbreak to understand the dynamics. In particular, we highlight that early data shows that the invading Omicron variant had a substantial growth advantage over the resident Delta variant. However, early data does not allow us to reliably infer other key epidemiological parameters—such as generation time and severity—which influence the expected peak hospital numbers. With more complete epidemic data from January 2022 are we able to capture the true scale of the epidemic in terms of both infections and hospital admissions, driven by different infection characteristics of Omicron compared to Delta and a substantial shift in estimated precautionary behaviour during December. This work highlights the challenges of real time forecasting, in a rapidly changing environment with limited information on the variant's epidemiological characteristics. Author summary: One of the key challenges for modellers at the start of any outbreak is to predict its likely scale and the associated burden on the health systems. This was the scenario in December of 2021, when the UK was faced with the invasion of the Omicron (BA.1) variant. For the Omicron BA.1 wave the peak in hospital admission was initially overestimated due to the simplifying assumption that Omicron would behave like a faster spreading version of the Delta variant, and that there would be little or no change to population behaviour during the wave. In this paper, we use a range of models and analyses to understand the challenges of correctly predicting the peak. Compared with forecasts made from the earliest data, the peak was reduced by the shorter infectious period of Omicron, its lower severity but mostly by changes in behaviour and testing over the Christmas holiday season. Historic projections made in December of 2021 that assumed a set of imposed mitigation measures were therefore more accurate than those that assumed no controls, although in practice the controls were largely self-imposed. This paper therefore highlights the challenges of forecasting the likely scale of outbreaks and the importance of understanding the population response to public health messaging. [ABSTRACT FROM AUTHOR]

In nature, most microbial populations have complex spatial structures that can affect their evolution. Evolutionary graph theory predicts that some spatial structures modelled by placing individuals on the nodes of a graph affect the probability that a mutant will fix. Evolution experiments are beginning to explicitly address the impact of graph structures on mutant fixation. However, the assumptions of evolutionary graph theory differ from the conditions of modern evolution experiments, making the comparison between theory and experiment challenging. Here, we aim to bridge this gap by using our new model of spatially structured populations. This model considers connected subpopulations that lie on the nodes of a graph, and allows asymmetric migrations. It can handle large populations, and explicitly models serial passage events with migrations, thus closely mimicking experimental conditions. We analyze recent experiments in light of this model. We suggest useful parameter regimes for future experiments, and we make quantitative predictions for these experiments. In particular, we propose experiments to directly test our recent prediction that the star graph with asymmetric migrations suppresses natural selection and can accelerate mutant fixation or extinction, compared to a well-mixed population. Author summary: Predicting how mutations spread through a population and eventually take over is important for understanding evolution. Complex spatial structures are ubiquitous in natural microbial populations, and can impact the fate of mutants. Theoretical models have been developed to describe this effect. They predict that some spatial structures have mutant fixation probabilities that differ from those of well-mixed populations. Experiments are beginning to probe these effects in the laboratory. However, there is a disconnect between models and experiments, because they consider different conditions. In this work, we connect them through a new model that closely matches experimental conditions. We analyze recent experiments and propose new ones that should allow testing the effects of complex population spatial structures on mutant fate. [ABSTRACT FROM AUTHOR]

Ligand binding to death receptors activates apoptosis in cancer cells. Stimulation of death receptors results in the formation of intracellular multiprotein platforms that either activate the apoptotic initiator Caspase-8 to trigger cell death, or signal through kinases to initiate inflammatory and cell survival signalling. Two of these platforms, the Death-Inducing Signalling Complex (DISC) and the RIPoptosome, also initiate necroptosis by building filamentous scaffolds that lead to the activation of mixed lineage kinase domain-like pseudokinase. To explain cell decision making downstream of death receptor activation, we developed a semi-stochastic model of DISC/RIPoptosome formation. The model is a hybrid of a direct Gillespie stochastic simulation algorithm for slow assembly of the RIPoptosome and a deterministic model of downstream caspase activation. The model explains how alterations in the level of death receptor-ligand complexes, their clustering properties and intrinsic molecular fluctuations in RIPoptosome assembly drive heterogeneous dynamics of Caspase-8 activation. The model highlights how kinetic proofreading leads to heterogeneous cell responses and results in fractional cell killing at low levels of receptor stimulation. It reveals that the noise in Caspase-8 activation—exclusively caused by the stochastic molecular assembly of the DISC/RIPoptosome platform—has a key function in extrinsic apoptotic stimuli recognition. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

Dynamic models in disease ecology have historically evaluated vaccination strategies under the assumption that they are implemented homogeneously in space and time. However, this approach fails to formally account for operational and logistical constraints inherent in the distribution of vaccination to the population at risk. Thus, feedback between the dynamic processes of vaccine distribution and transmission might be overlooked. Here, we present a spatially explicit, stochastic Susceptible-Infected-Recovered-Vaccinated model that highlights the density-dependence and spatial constraints of various diffusive strategies of vaccination during an outbreak. The model integrates an agent-based process of disease spread with a partial differential process of vaccination deployment. We characterize the vaccination response in terms of a diffusion rate that describes the distribution of vaccination to the population at risk from a central location. This generates an explicit trade-off between slow diffusion, which concentrates effort near the central location, and fast diffusion, which spreads a fixed vaccination effort thinly over a large area. We use stochastic simulation to identify the optimum vaccination diffusion rate as a function of population density, interaction scale, transmissibility, and vaccine intensity. Our results show that, conditional on a timely response, the optimal strategy for minimizing outbreak size is to distribute vaccination resource at an intermediate rate: fast enough to outpace the epidemic, but slow enough to achieve local herd immunity. If the response is delayed, however, the optimal strategy for minimizing outbreak size changes to a rapidly diffusive distribution of vaccination effort. The latter may also result in significantly larger outbreaks, thus suggesting a benefit of allocating resources to timely outbreak detection and response. [ABSTRACT FROM AUTHOR]

Reducing spillover of zoonotic pathogens is an appealing approach to preventing human disease and minimizing the risk of future epidemics and pandemics. Although the immediate human health benefit of reducing spillover is clear, over time, spillover reduction could lead to counterintuitive negative consequences for human health. Here, we use mathematical models and computer simulations to explore the conditions under which unanticipated consequences of spillover reduction can occur in systems where the severity of disease increases with age at infection. Our results demonstrate that, because the average age at infection increases as spillover is reduced, programs that reduce spillover can actually increase population-level disease burden if the clinical severity of infection increases sufficiently rapidly with age. If, however, immunity wanes over time and reinfection is possible, our results reveal that negative health impacts of spillover reduction become substantially less likely. When our model is parameterized using published data on Lassa virus in West Africa, it predicts that negative health outcomes are possible, but likely to be restricted to a small subset of populations where spillover is unusually intense. Together, our results suggest that adverse consequences of spillover reduction programs are unlikely but that the public health gains observed immediately after spillover reduction may fade over time as the age structure of immunity gradually re-equilibrates to a reduced force of infection. Author summary: Many pathogens, such as rabies, coronaviruses, and hantaviruses primarily circulate within wild animals but can infect humans when the opportunity arises. This pervasive challenge to public health has motivated the development of new methods designed to reduce the frequency of these spillover events. Although reducing spillover infection of humans appears to be an obvious win for public health, it is conceivable that altering historical patterns of spillover could change the age structure of human immunity in a way that undermines human health. Using mathematical and computational models we evaluate the conditions required for these counterintuitive impacts to occur. Our analyses demonstrate that reducing spillover from wild animals will generally improve public health and that negative outcomes can occur in only rare and unusual circumstances. Although negative impacts of spillover reduction are likely to be rare, our results show that the public health benefits of spillover reduction may fade over time unless a near total elimination of spillover can be achieved. [ABSTRACT FROM AUTHOR]

Hepatic stellate cells (HSC) are the source of extracellular matrix (ECM) whose overproduction leads to fibrosis, a condition that impairs liver functions in chronic liver diseases. Understanding the dynamics of HSCs will provide insights needed to develop new therapeutic approaches. Few models of hepatic fibrosis have been proposed, and none of them include the heterogeneity of HSC phenotypes recently highlighted by single-cell RNA sequencing analyses. Here, we developed rule-based models to study HSC dynamics during fibrosis progression and reversion. We used the Kappa graph rewriting language, for which we used tokens and counters to overcome temporal explosion. HSCs are modeled as agents that present seven physiological cellular states and that interact with (TGFβ1) molecules which regulate HSC activation and the secretion of type I collagen, the main component of the ECM. Simulation studies revealed the critical role of the HSC inactivation process during fibrosis progression and reversion. While inactivation allows elimination of activated HSCs during reversion steps, reactivation loops of inactivated HSCs (iHSCs) are required to sustain fibrosis. Furthermore, we demonstrated the model's sensitivity to (TGFβ1) parameters, suggesting its adaptability to a variety of pathophysiological conditions for which levels of (TGFβ1) production associated with the inflammatory response differ. Using new experimental data from a mouse model of CCl4-induced liver fibrosis, we validated the predicted ECM dynamics. Our model also predicts the accumulation of iHSCs during chronic liver disease. By analyzing RNA sequencing data from patients with non-alcoholic steatohepatitis (NASH) associated with liver fibrosis, we confirmed this accumulation, identifying iHSCs as novel markers of fibrosis progression. Overall, our study provides the first model of HSC dynamics in chronic liver disease that can be used to explore the regulatory role of iHSCs in liver homeostasis. Moreover, our model can also be generalized to fibroblasts during repair and fibrosis in other tissues. Author summary: Chronic liver diseases (CLDs) are associated with the development of fibrosis which is characterized by an abnormal deposition of extracellular matrix (ECM) leading to severe liver dysfunction. Hepatic stellate cells (HSCs) are key players in liver fibrosis driving ECM remodeling. However numerous biological processes are involved including HSC activation, proliferation, differentiation and inactivation and novel computational modeling is necessary to integrate such complex dynamics. Here, we used the Kappa graph rewriting language to develop the first rule-based model describing the HSCs dynamics during liver fibrosis and its reversion. Simulation analyses enabled us to demonstrate the critical role of the HSC inactivation loop in the development of liver fibrosis, and to identify inactivated HSCs as potential new markers of fibrosis progression. [ABSTRACT FROM AUTHOR]

Cells switch genes ON or OFF by altering the state of chromatin via histone modifications at specific regulatory locations along the chromatin polymer. These gene regulation processes are carried out by a network of reactions in which the histone marks spread to neighboring regions with the help of enzymes. In the literature, this spreading has been studied as a purely kinetic, non-diffusive process considering the interactions between neighboring nucleosomes. In this work, we go beyond this framework and study the spreading of modifications using a reaction-diffusion (RD) model accounting for the diffusion of the constituents. We quantitatively segregate the modification profiles generated from kinetic and RD models. The diffusion and degradation of enzymes set a natural length scale for limiting the domain size of modification spreading, and the resulting enzyme limitation is inherent in our model. We also demonstrate the emergence of confined modification domains without the explicit requirement of a nucleation site. We explore polymer compaction effects on spreading and show that single-cell domains may differ from averaged profiles. We find that the modification profiles from our model are comparable with existing H3K9me3 data of S. pombe. Author summary: DNA is folded into chromatin by wrapping it around histone protein octamers, forming a string of nucleosome beads. The chemical modifications on histone proteins along the one-dimensional string form a significant part of the epigenetic code. This code determines whether a chromatin region is likely to be active or repressed. How the chemical modification patterns arise, and their stability along the histone proteins is a question of great interest. While some players in this process are known, how the physical process of diffusion and reaction of the enzymes collectively dictate the pattern formation is unknown. In this work, we simulate histone modification spreading accounting for explicit diffusion and reaction of the essential constituents. We show that the diffusion and degradation processes of the constituents limit the size of the modification domain. We also show how the availability of enzymes and the compaction of the chromatin affect the spreading of the histone modification. [ABSTRACT FROM AUTHOR]

Amortized simulation-based neural posterior estimation provides a novel machine learning based approach for solving parameter estimation problems. It has been shown to be computationally efficient and able to handle complex models and data sets. Yet, the available approach cannot handle the in experimental studies ubiquitous case of missing data, and might provide incorrect posterior estimates. In this work, we discuss various ways of encoding missing data and integrate them into the training and inference process. We implement the approaches in the BayesFlow methodology, an amortized estimation framework based on invertible neural networks, and evaluate their performance on multiple test problems. We find that an approach in which the data vector is augmented with binary indicators of presence or absence of values performs the most robustly. Indeed, it improved the performance also for the simpler problem of data sets with variable length. Accordingly, we demonstrate that amortized simulation-based inference approaches are applicable even with missing data, and we provide a guideline for their handling, which is relevant for a broad spectrum of applications. Author summary: In biomedical research, mechanistic models describe dynamic processes, yet inferring their underlying parameters can often be challenging. Bayesian statistics provides an established framework for this by integrating prior knowledge with observed data, and naturally enables uncertainty quantification as a distribution of parameter values is returned. However, classical case-based methods for Bayesian inference can be computationally expensive, particularly when the same model needs to be fitted to different data sets. Recently, deep-learning-based approaches have been developed to streamline the inference procedure, allowing the upfront training cost to amortize when applied to multiple data sets. In this manuscript, we explore approaches to extend the setup to data sets with missing data. In summary, an encoding scheme which exploits data augmentation with binary indicators of presence or absence performs the most robustly across different test problems. [ABSTRACT FROM AUTHOR]

Contagion processes, representing the spread of infectious diseases, information, or social behaviors, are often schematized as taking place on networks, which encode for instance the interactions between individuals. The impact of the network structure on spreading process has been widely investigated, but not the reverse question: do different processes unfolding on a given network lead to different infection patterns? How do the infection patterns depend on a model's parameters or on the nature of the contagion processes? Here we address this issue by investigating the infection patterns for a variety of models. In simple contagion processes, where contagion events involve one connection at a time, we find that the infection patterns are extremely robust across models and parameters. In complex contagion models instead, in which multiple interactions are needed for a contagion event, non-trivial dependencies on models parameters emerge, as the infection pattern depends on the interplay between pairwise and group contagions. In models involving threshold mechanisms moreover, slight parameter changes can significantly impact the spreading paths. Our results show that it is possible to study crucial features of a spread from schematized models, and inform us on the variations between spreading patterns in processes of different nature. Author summary: Contagion processes, representing the spread of infectious diseases, information, or social behaviors, are often schematized as taking place on networks, which encode for instance the interactions between individuals. We here observe how the network is explored by the contagion process, i.e. which links are used for contagions and how frequently. The resulting infection pattern depends on the chosen infection model but surprisingly not all the parameters and models features play a role in the infection pattern. We discover for instance that in simple contagion processes, where contagion events involve one connection at a time, the infection patterns are extremely robust across models and parameters. This has consequences in the role of models in decision-making, as it implies that numerical simulations of simple contagion processes using simplified settings can bring important insights even in the case of a new emerging disease whose properties are not yet well known. In complex contagion models instead, in which multiple interactions are needed for a contagion event, non-trivial dependencies on model parameters emerge and infection pattern cannot be confused with those observed for simple contagion. [ABSTRACT FROM AUTHOR]

Cellular signaling, crucial for biological processes like immune response and homeostasis, relies on specificity and fidelity in signal transduction to accurately respond to stimuli amidst biological noise. Kinetic proofreading (KPR) is a key mechanism enhancing signaling specificity through time-delayed steps, although its effectiveness is debated due to intrinsic noise potentially reducing signal fidelity. In this study, we reformulate the theory of kinetic proofreading (KPR) by convolving multiple intermediate states into a single state and then define an overall "processing" time required to traverse these states. This simplification allows us to succinctly describe kinetic proofreading in terms of a single waiting time parameter, facilitating a more direct evaluation and comparison of KPR performance across different biological contexts such as DNA replication and T cell receptor (TCR) signaling. We find that loss of fidelity for longer proofreading steps relies on the specific strategy of information extraction and show that in the first-passage time (FPT) discrimination strategy, longer proofreading steps can exponentially improve the accuracy of KPR at the cost of speed. Thus, KPR can still be an effective discrimination mechanism in the high noise regime. However, in a product concentration-based discrimination strategy, longer proofreading steps do not necessarily lead to an increase in performance. However, by introducing activation thresholds on product concentrations, can we decompose the product-based strategy into a series of FPT-based strategies to better resolve the subtleties of KPR-mediated product discrimination. Our findings underscore the importance of understanding KPR in the context of how information is extracted and processed in the cell. Author summary: Kinetic proofreading (KPR) is mechanism often employed by cells to enhance specificity of ligand-receptor. However, the performance of kinetic proofreading may be hampered by noise and a low signal-to-noise ratio. By consolidating multiple kinetic proofreading steps into a single state and assigning an associated waiting, or "processing time," we developed an analytic approach to quantify the performance of KPR in different biological contexts. Despite a trade-off between speed and accuracy inherent to a first-passage time KPR strategy, we show that a signaling molecule-based discrimination strategy can enhance the performance benefits of KPR. We further decompose the product-based discrimination strategy into a set of first-passage times to different thresholds of signaling molecules produced. Through this decomposition, we find that a threshold that adjusts dynamically throughout the recognition process depends on the duration of the process. We propose that this more nuanced product-based KPR-mediated recognition process can be realized biologically. The precise structural basis for a dynamic threshold merits further experimental exploration, as it may hold significant implications for understanding biological mechanisms of information transmission at a molecular level. [ABSTRACT FROM AUTHOR]

The Michaelis–Menten (MM) rate law has been a fundamental tool in describing enzyme-catalyzed reactions for over a century. When substrates and enzymes are homogeneously distributed, the validity of the MM rate law can be easily assessed based on relative concentrations: the substrate is in large excess over the enzyme-substrate complex. However, the applicability of this conventional criterion remains unclear when species exhibit spatial heterogeneity, a prevailing scenario in biological systems. Here, we explore the MM rate law's applicability under spatial heterogeneity by using partial differential equations. In this study, molecules diffuse very slowly, allowing them to locally reach quasi-steady states at each time point. We find that the conventional criterion for the validity of the MM rate law cannot be readily extended to heterogeneous environments solely through spatial averages of molecular concentrations. That is, even when the conventional criterion for the spatial averages is satisfied, the MM rate law fails to capture the enzyme catalytic rate under spatial heterogeneity. In contrast, a slightly modified form of the MM rate law, based on the total quasi-steady state approximation (tQSSA), is accurate. Specifically, the tQSSA-based modified form, but not the original MM rate law, accurately predicts the drug clearance via cytochrome P450 enzymes and the ultrasensitive phosphorylation in heterogeneous environments. Our findings shed light on how to simplify spatiotemporal models for enzyme-catalyzed reactions in the right context, ensuring accurate conclusions and avoiding misinterpretations in in silico simulations. Author summary: For over a century, scientists have relied on the simple Michaelis–Menten (MM) rate law to explain how enzymes function. The conventional criterion for using the MM rate law has been derived in homogeneous environments where species are evenly dispersed (e.g., in vitro experiments). Here, we found that the conventional criterion does not work in intracellular environments where species are heterogeneously distributed and diffuse slowly. However, we find that a slightly modified formula of the MM rate law, based on total quasi-steady state approximation (tQSSA), is accurate even in heterogeneous environments. In particular, this modified formula, unlike the MM rate law, accurately predicts the rate of drug metabolism and ultrasensitive phosphorylation when species are not evenly distributed. Our results provide insight into how to use simplified models for describing enzyme functions in such non-evenly distributed environments. [ABSTRACT FROM AUTHOR]

Pre-exposure prophylaxis (PrEP) is an important pillar to prevent HIV transmission. Because of experimental and clinical shortcomings, mathematical models that integrate pharmacological, viral- and host factors are frequently used to quantify clinical efficacy of PrEP. Stochastic simulations of these models provides sample statistics from which the clinical efficacy is approximated. However, many stochastic simulations are needed to reduce the associated sampling error. To remedy the shortcomings of stochastic simulation, we developed a numerical method that allows predicting the efficacy of arbitrary prophylactic regimen directly from a viral dynamics model, without sampling. We apply the method to various hypothetical dolutegravir (DTG) prophylaxis scenarios. The approach is verified against state-of-the-art stochastic simulation. While the method is more accurate than stochastic simulation, it is superior in terms of computational performance. For example, a continuous 6-month prophylactic profile is computed within a few seconds on a laptop computer. The method's computational performance, therefore, substantially expands the horizon of feasible analysis in the context of PrEP, and possibly other applications.

Exposure to environmental stressors, including certain antibiotics, induces stress responses in bacteria. Some of these responses increase mutagenesis and thus potentially accelerate resistance evolution. Many studies report increased mutation rates under stress, often using the standard experimental approach of fluctuation assays. However, single-cell studies have revealed that many stress responses are heterogeneously expressed in bacterial populations, which existing estimation methods have not yet addressed. We develop a population dynamic model that considers heterogeneous stress responses (subpopulations of cells with the response off or on) that impact both mutation rate and cell division rate, inspired by the DNA-damage response in Escherichia coli (SOS response). We derive the mutant count distribution arising in fluctuation assays under this model and then implement maximum likelihood estimation of the mutation-rate increase specifically associated with the expression of the stress response. Using simulated mutant count data, we show that our inference method allows for accurate and precise estimation of the mutation-rate increase, provided that this increase is sufficiently large and the induction of the response also reduces the division rate. Moreover, we find that in many cases, either heterogeneity in stress responses or mutant fitness costs could explain similar patterns in fluctuation assay data, suggesting that separate experiments would be required to identify the true underlying process. In cases where stress responses and mutation rates are heterogeneous, current methods still correctly infer the effective increase in population mean mutation rate, but we provide a novel method to infer distinct stress-induced mutation rates, which could be important for parameterising evolutionary models. Author summary: How does environmental stress, especially from antibiotics, affect mutation rates in bacteria? This question has often been examined by estimating mutation rates using fluctuation assays, an experiment dating back to Luria and Delbrück in the 1940s. In this study, we consider variation in stress responses within bacterial populations, as revealed by recent single-cell studies, which is neglected in currently available mutation-rate estimation methods. Our approach involves a population dynamic model inspired by the DNA-damage response in E. coli (SOS response). It accounts for a subpopulation with high expression of the stress response, which increases the mutation rate and decreases the division rate of a cell. We use computer simulations to generate synthetic fluctuation assay data. Notably, we find that over a wide range of scenarios, existing models and our heterogeneous-response model cannot be distinguished using fluctuation assay data alone. This emphasises the need for separate experiments to uncover the true underlying processes. Nevertheless, when stress responses are known to be heterogeneous, our study offers a novel method for accurately estimating mutation rates specifically associated with the high expression of the stress response. Uncovering the heterogeneity in stress-induced mutation rates could be important for predicting the evolution of antibiotic resistance. [ABSTRACT FROM AUTHOR]

The self-organization of cells relies on the profound complexity of protein-protein interactions. Challenges in directly observing these events have hindered progress toward understanding their diverse behaviors. One notable example is the interaction between molecular motors and cytoskeletal systems that combine to perform a variety of cellular functions. In this work, we leverage theory and experiments to identify and quantify the rate-limiting mechanism of the initial association between a cargo-bound kinesin motor and a microtubule track. Recent advances in optical tweezers provide binding times for several lengths of kinesin motors trapped at varying distances from a microtubule, empowering the investigation of competing models. We first explore a diffusion-limited model of binding. Through Brownian dynamics simulations and simulation-based inference, we find this simple diffusion model fails to explain the experimental binding times, but an extended model that accounts for the ADP state of the molecular motor agrees closely with the data, even under the scrutiny of penalizing for additional model complexity. We provide quantification of both kinetic rates and biophysical parameters underlying the proposed binding process. Our model suggests that a typical binding event is limited by ADP state rather than physical search. Lastly, we predict how these association rates can be modulated in distinct ways through variation of environmental concentrations and physical properties. Author summary: Cytoskeletal-motor assemblies self-organize to achieve cellular functions ranging from delivering intracellular cargoes to generating forces in mitosis. Advancements in single-molecule experiments have revealed immense detail about motor detachment and stepping, but relatively little regarding attachment. Newly available binding times for individual kinesin motors allow for the evaluation of mechanistic models of the process. We find that a model limited by both diffusive search and a weak-to-strong transition from ADP release best explains the data. The coupled chemo-mechanical nature of this interaction can be modulated more richly than either separately, possibly explaining the diversity and regulation observed in cells. More broadly, our study provides a timely vignette in leveraging computations with experiments to understand the mosaic of mechanisms underlying protein-protein interactions. [ABSTRACT FROM AUTHOR]

In this work, we quantitatively compare computer simulations and existing cell tracking data of P. aeruginosa surface motility in order to analyse the underlying motility mechanism. We present a three dimensional twitching motility model, that simulates the extension, retraction and surface association of individual Type IV Pili (TFP), and is informed by recent experimental observations of TFP. Sensitivity analysis is implemented to minimise the number of model parameters, and quantitative estimates for the remaining parameters are inferred from tracking data by approximate Bayesian computation. We argue that the motility mechanism is highly sensitive to experimental conditions. We predict a TFP retraction speed for the tracking data we study that is in a good agreement with experimental results obtained under very similar conditions. Furthermore, we examine whether estimates for biologically important parameters, whose direct experimental determination is challenging, can be inferred directly from tracking data. One example is the width of the distribution of TFP on the bacteria body. We predict that the TFP are broadly distributed over the bacteria pole in both walking and crawling motility types. Moreover, we identified specific configurations of TFP that lead to transitions between walking and crawling states. Author summary: Twitching is a type of bacterial surface motion facilitated by molecular motor-driven micron-scale filaments known as type IV pili (TFP). The resulting motion, which is characterised by erratic behaviour at short timescales, is important for bacterial surface invasion, biofilm formation etc. The link between microscopic mechanisms on the level of single filaments and macroscopic properties of twitching trajectories is generally difficult to make, since direct time-resolved imaging is of TFP in experiments is challenging. We propose a data-driven model for surface motility of the bacteria P. Aeruginosa that bridges these scales: we resolve single-TFP dynamics and connect it to the macroscopic properties of twitching trajectories. We quantitatively compare our simulations to experimental tracking data, predict the previously unresolved distribution of TFP on bacterial membrane, and discuss how it affects the transitions between the "walking" and "crawling" modes of motion. We also extract the speed of TFP retraction corresponding to the specific experimental conditions, which is lower than reported in several other experiments, but in perfect agreement with measurements performed at conditions of our tracking data. Our work thus elucidates microscopic mechanism of twitching, and it represents a systematic data-driven approach to quantitative modelling of complex biological phenomena. [ABSTRACT FROM AUTHOR]

In experiments, the distributions of mRNA or protein numbers in single cells are often fitted to the random telegraph model which includes synthesis and decay of mRNA or protein, and switching of the gene between active and inactive states. While commonly used, this model does not describe how fluctuations are influenced by crucial biological mechanisms such as feedback regulation, non-exponential gene inactivation durations, and multiple gene activation pathways. Here we investigate the dynamical properties of four relatively complex gene expression models by fitting their steady-state mRNA or protein number distributions to the simple telegraph model. We show that despite the underlying complex biological mechanisms, the telegraph model with three effective parameters can accurately capture the steady-state gene product distributions, as well as the conditional distributions in the active gene state, of the complex models. Some effective parameters are reliable and can reflect realistic dynamic behaviors of the complex models, while others may deviate significantly from their real values in the complex models. The effective parameters can also be applied to characterize the capability for a complex model to exhibit multimodality. Using additional information such as single-cell data at multiple time points, we provide an effective method of distinguishing the complex models from the telegraph model. Furthermore, using measurements under varying experimental conditions, we show that fitting the mRNA or protein number distributions to the telegraph model may even reveal the underlying gene regulation mechanisms of the complex models. The effectiveness of these methods is confirmed by analysis of single-cell data for E. coli and mammalian cells. All these results are robust with respect to cooperative transcriptional regulation and extrinsic noise. In particular, we find that faster relaxation speed to the steady state results in more precise parameter inference under large extrinsic noise. Author summary: Over the past decade, significant progress has been made in the theory and experiments of single-cell stochastic gene expression dynamics. The most well studied and widely used stochastic gene expression model is the two-state telegraph model. However, the conventional telegraph model is too simple and limited in its predictive power because it lacks a description of some important biological mechanisms, such as feedback regulation, multiple gene activation steps, and multiple gene activation pathways. This raises a important question: what can we learn when fitting a complex gene expression model to a simple telegraph model? In this paper, we investigate four complex gene expression models by fitting their steady-state mRNA or protein number distributions to the telegraph model and then obtain estimates of the effective parameters. We show that while the estimated values of the parameters in the "artificial" telegraph model are not always accurate, they are still sometimes reliable and can also reveal important dynamical properties of the complex models such as the ability for a complex model to produce bimodality. Moreover, we provide an effective method of distinguishing the complex models from the telegraph model by using additional information such as gene expression data at multiple time points. Finally, we show that fitting the mRNA or protein number distributions to the telegraph model may even reveal the underlying gene regulation mechanism of a complex model by using measurements under varying experimental conditions. The effectiveness of these methods is well confirmed by analysis of single-cell gene expression data for E. coli and mammalian cells. [ABSTRACT FROM AUTHOR]

Biochemical signaling pathways in living cells are often highly organized into spatially segregated volumes, membranes, scaffolds, subcellular compartments, and organelles comprising small numbers of interacting molecules. At this level of granularity stochastic behavior dominates, well-mixed continuum approximations based on concentrations break down and a particle-based approach is more accurate and more efficient. We describe and validate a new version of the open-source MCell simulation program (MCell4), which supports generalized 3D Monte Carlo modeling of diffusion and chemical reaction of discrete molecules and macromolecular complexes in solution, on surfaces representing membranes, and combinations thereof. The main improvements in MCell4 compared to the previous versions, MCell3 and MCell3-R, include a Python interface and native BioNetGen reaction language (BNGL) support. MCell4's Python interface opens up completely new possibilities for interfacing with external simulators to allow creation of sophisticated event-driven multiscale/multiphysics simulations. The native BNGL support, implemented through a new open-source library libBNG (also introduced in this paper), provides the capability to run a given BNGL model spatially resolved in MCell4 and, with appropriate simplifying assumptions, also in the BioNetGen simulation environment, greatly accelerating and simplifying model validation and comparison. [ABSTRACT FROM AUTHOR]

In the midst of an outbreak or sustained epidemic, reliable prediction of transmission risks and patterns of spread is critical to inform public health programs. Projections of transmission growth or decline among specific risk groups can aid in optimizing interventions, particularly when resources are limited. Phylogenetic trees have been widely used in the detection of transmission chains and high-risk populations. Moreover, tree topology and the incorporation of population parameters (phylodynamics) can be useful in reconstructing the evolutionary dynamics of an epidemic across space and time among individuals. We now demonstrate the utility of phylodynamic trees for transmission modeling and forecasting, developing a phylogeny-based deep learning system, referred to as DeepDynaForecast. Our approach leverages a primal-dual graph learning structure with shortcut multi-layer aggregation, which is suited for the early identification and prediction of transmission dynamics in emerging high-risk groups. We demonstrate the accuracy of DeepDynaForecast using simulated outbreak data and the utility of the learned model using empirical, large-scale data from the human immunodeficiency virus epidemic in Florida between 2012 and 2020. Our framework is available as open-source software (MIT license) at github.com/lab-smile/DeepDynaForcast. Author summary: During an outbreak or sustained epidemic, accurate prediction of patterns in transmission risk can reliably inform public health strategies. Projections indicating growth or decline of transmission for specific risk groups can significantly enhance the optimization of interventions, especially when resources are limited. To address this, we present DeepDynaForecast, a cutting-edge deep learning algorithm designed for forecasting pathogen transmission dynamics. Uniquely, DeepDynaForecast was trained on in-depth simulation data, classifying samples according to their dynamics (growth, static, or decline) with accuracy of 91.6%. We evaluated the model's performance and application using simulated outbreak data and empirical, large-scale data from the HIV epidemic in Florida between 2012 and 2020. We conclude DeepDynaForecast represents a significant advancement in genomics-mediated pathogen transmission characterization and has the potential to catalyze new research directions within virology, molecular biology, and public health. [ABSTRACT FROM AUTHOR]