Systematic comparison of modeling fidelity levels and parameter inference settings applied to negative feedback gene regulation.

Quantitative stochastic models of gene regulatory networks are important tools for studying cellular regulation. Such models can be formulated at many different levels of fidelity. A practical challenge is to determine what model fidelity to use in order to get accurate and representative results. The choice is important, because models of successively higher fidelity come at a rapidly increasing computational cost. In some situations, the level of detail is clearly motivated by the question under study. In many situations however, many model options could qualitatively agree with available data, depending on the amount of data and the nature of the observations. Here, an important distinction is whether we are interested in inferring the true (but unknown) physical parameters of the model or if it is sufficient to be able to capture and explain available data. The situation becomes complicated from a computational perspective because inference needs to be approximate. Most often it is based on likelihood-free Approximate Bayesian Computation (ABC) and here determining which summary statistics to use, as well as how much data is needed to reach the desired level of accuracy, are difficult tasks. Ultimately, all of these aspects-the model fidelity, the available data, and the numerical choices for inference-interplay in a complex manner. In this paper we develop a computational pipeline designed to systematically evaluate inference accuracy for a wide range of true known parameters. We then use it to explore inference settings for negative feedback gene regulation. In particular, we compare a detailed spatial stochastic model, a coarse-grained compartment-based multiscale model, and the standard well-mixed model, across several data-scenarios and for multiple numerical options for parameter inference. Practically speaking, this pipeline can be used as a preliminary step to guide modelers prior to gathering experimental data. By training Gaussian processes to approximate the distance function values, we are able to substantially reduce the computational cost of running the pipeline.

Efficient Bayesian inference for stochastic agent-based models.

The modelling of many real-world problems relies on computationally heavy simulations of randomly interacting individuals or agents. However, the values of the parameters that underlie the interactions between agents are typically poorly known, and hence they need to be inferred from macroscopic observations of the system. Since statistical inference rests on repeated simulations to sample the parameter space, the high computational expense of these simulations can become a stumbling block. In this paper, we compare two ways to mitigate this issue in a Bayesian setting through the use of machine learning methods: One approach is to construct lightweight surrogate models to substitute the simulations used in inference. Alternatively, one might altogether circumvent the need for Bayesian sampling schemes and directly estimate the posterior distribution. We focus on stochastic simulations that track autonomous agents and present two case studies: tumour growths and the spread of infectious diseases. We demonstrate that good accuracy in inference can be achieved with a relatively small number of simulations, making our machine learning approaches orders of magnitude faster than classical simulation-based methods that rely on sampling the parameter space. However, we find that while some methods generally produce more robust results than others, no algorithm offers a one-size-fits-all solution when attempting to infer model parameters from observations. Instead, one must choose the inference technique with the specific real-world application in mind. The stochastic nature of the considered real-world phenomena poses an additional challenge that can become insurmountable for some approaches. Overall, we find machine learning approaches that create direct inference machines to be promising for real-world applications. We present our findings as general guidelines for modelling practitioners.

Benchmarking imputation methods for network inference using a novel method of synthetic scRNA-seq data generation.

BACKGROUND: Single cell RNA-sequencing (scRNA-seq) has very rapidly become the new workhorse of modern biology providing an unprecedented global view on cellular diversity and heterogeneity. In particular, the structure of gene-gene expression correlation contains information on the underlying gene regulatory networks. However, interpretation of scRNA-seq data is challenging due to specific experimental error and biases that are unique to this kind of data including drop-out (or technical zeros). METHODS: To deal with this problem several methods for imputation of zeros for scRNA-seq have been developed. However, it is not clear how these processing steps affect inference of genetic networks from single cell data. Here, we introduce Biomodelling.jl, a tool for generation of synthetic scRNA-seq data using multiscale modelling of stochastic gene regulatory networks in growing and dividing cells. RESULTS: Our tool produces realistic transcription data with a known ground truth network topology that can be used to benchmark different approaches for gene regulatory network inference. Using this tool we investigate the impact of different imputation methods on the performance of several network inference algorithms. CONCLUSIONS: Biomodelling.jl provides a versatile and useful tool for future development and benchmarking of network inference approaches using scRNA-seq data.

The influence of methylation status on a stochastic model of MGMT dynamics in glioblastoma: Phenotypic selection can occur with and without a downshift in promoter methylation status.

Glioblastoma originates in the brain and is one of the most aggressive cancer types. Glioblastoma represents 15% of all brain tumours, with a median survival of 15 months. Although the current standard of care for such a tumour (the Stupp protocol) has shown positive results for the prognosis of patients, O-6-methylguanine-DNA methyltransferase (MGMT) driven drug resistance has been an issue of increasing concern and hence requires innovative approaches. In addition to the well established drug resistance factors such as tumour location and blood brain barriers, it is also important to understand how the genetic and epigenetic dynamics of the glioblastoma cells can play a role. One important aspect of this is the study of methylation status of MGMT following administration of temozolomide. In this paper, we extend our previously published model that simulated MGMT expression in glioblastoma cells to incorporate the promoter methylation status of MGMT. This methylation status has clinical significance and is used as a marker for patient outcomes. Using this model, we investigate the causative relationship between temozolomide treatment and the methylation status of the MGMT promoter in a population of cells. In addition by constraining the model to relevant biological data using Approximate Bayesian Computation, we were able to identify parameter regimes that yield different possible modes of resistances, namely, phenotypic selection of MGMT, a downshift in the methylation status of the MGMT promoter or both simultaneously. We analysed each of the parameter sets associated with the different modes of resistance, presenting representative solutions as well as discovering some similarities between them as well as unique requirements for each of them. Finally, we used them to devise optimal strategies for inhibiting MGMT expression with the aim of minimising live glioblastoma cell numbers.

Heterogeneous responses to low level death receptor activation are explained by random molecular assembly of the Caspase-8 activation platform.

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.

The influence of nuclear compartmentalisation on stochastic dynamics of self-repressing gene expression.

Gene expression is an inherently noisy process. This noise is generally thought to be deleterious as precise internal regulation of biochemical reactions is essential for cell growth and survival. Self-repression of gene expression, which is the simplest form of a negative feedback loop, is commonly believed to be employed by cellular systems to decrease the stochastic fluctuations in gene expression. When there is some delay in autoregulation, it is also believed that this system can generate oscillations. In eukaryotic cells, mRNAs that are synthesised in the nucleus must be exported to the cytoplasm to function in protein synthesis, whereas proteins must be transported into the nucleus from the cytoplasm to regulate the expression levels of genes. Nuclear transport thus plays a critical role in eukaryotic gene expression and regulation. Some recent studies have suggested that nuclear retention of mRNAs can control noise in mRNA expression. However, the effect of nuclear transport on protein noise and its interplay with negative feedback regulation is not completely understood. In this paper, we systematically compare four different simple models of gene expression. By using simulations and applying the linear noise approximation to the corresponding chemical master equations, we investigate the influence of nuclear import and export on noise in gene expression in a negative autoregulatory feedback loop. We first present results consistent with the literature, i.e., that negative feedback can effectively buffer the variability in protein levels, and nuclear retention can decrease mRNA noise levels. Interestingly we find that when negative feedback is combined with nuclear retention, an amplification in gene expression noise can be observed and is dependant on nuclear translocation rates. Finally, we investigate the effect of nuclear compartmentalisation on the ability of self-repressing genes to exhibit stochastic oscillatory dynamics.

Stochastic reaction-diffusion algorithms for macromolecular crowding.

Compartment-based (lattice-based) reaction-diffusion algorithms are often used for studying complex stochastic spatio-temporal processes inside cells. In this paper the influence of macromolecular crowding on stochastic reaction-diffusion simulations is investigated. Reaction-diffusion processes are considered on two different kinds of compartmental lattice, a cubic lattice and a hexagonal close packed lattice, and solved using two different algorithms, the stochastic simulation algorithm and the spatiocyte algorithm (Arjunan and Tomita 2010 Syst. Synth. Biol. 4, 35-53). Obstacles (modelling macromolecular crowding) are shown to have substantial effects on the mean squared displacement and average number of molecules in the domain but the nature of these effects is dependent on the choice of lattice, with the cubic lattice being more susceptible to the effects of the obstacles. Finally, improvements for both algorithms are presented.

A mathematical model of pre-diagnostic glioma growth.

Due to their location, the malignant gliomas of the brain in humans are very difficult to treat in advanced stages. Blood-based biomarkers for glioma are needed for more accurate evaluation of treatment response as well as early diagnosis. However, biomarker research in primary brain tumors is challenging given their relative rarity and genetic diversity. It is further complicated by variations in the permeability of the blood brain barrier that affects the amount of marker released into the bloodstream. Inspired by recent temporal data indicating a possible decrease in serum glucose levels in patients with gliomas yet to be diagnosed, we present an ordinary differential equation model to capture early stage glioma growth. The model contains glioma-glucose-immune interactions and poses a potential mechanism by which this glucose drop can be explained. We present numerical simulations, parameter sensitivity analysis, linear stability analysis and a numerical experiment whereby we show how a dormant glioma can become malignant.

The role of dimerisation and nuclear transport in the Hes1 gene regulatory network.

Hes1 is a member of the family of basic helix-loop-helix transcription factors and the Hes1 gene regulatory network (GRN) may be described as the canonical example of transcriptional control in eukaryotic cells, since it involves only the Hes1 protein and its own mRNA. Recently, the Hes1 protein has been established as an excellent target for an anti-cancer drug treatment, with the design of a small molecule Hes1 dimerisation inhibitor representing a promising if challenging approach to therapy. In this paper, we extend a previous spatial stochastic model of the Hes1 GRN to include nuclear transport and dimerisation of Hes1 monomers. Initially, we assume that dimerisation occurs only in the cytoplasm, with only dimers being imported into the nucleus. Stochastic simulations of this novel model using the URDME software show that oscillatory dynamics in agreement with experimental studies are retained. Furthermore, we find that our model is robust to changes in the nuclear transport and dimerisation parameters. However, since the precise dynamics of the nuclear import of Hes1 and the localisation of the dimerisation reaction are not known, we consider a second modelling scenario in which we allow for both Hes1 monomers and dimers to be imported into the nucleus, and we allow dimerisation of Hes1 to occur everywhere in the cell. Once again, computational solutions of this second model produce oscillatory dynamics in agreement with experimental studies. We also explore sensitivity of the numerical solutions to nuclear transport and dimerisation parameters. Finally, we compare and contrast the two different modelling scenarios using numerical experiments that simulate dimer disruption, and suggest a biological experiment that could distinguish which model more faithfully captures the Hes1 GRN.

Data-driven spatio-temporal modelling of glioblastoma.

Mathematical oncology provides unique and invaluable insights into tumour growth on both the microscopic and macroscopic levels. This review presents state-of-the-art modelling techniques and focuses on their role in understanding glioblastoma, a malignant form of brain cancer. For each approach, we summarize the scope, drawbacks and assets. We highlight the potential clinical applications of each modelling technique and discuss the connections between the mathematical models and the molecular and imaging data used to inform them. By doing so, we aim to prime cancer researchers with current and emerging computational tools for understanding tumour progression. By providing an in-depth picture of the different modelling techniques, we also aim to assist researchers who seek to build and develop their own models and the associated inference frameworks. Our article thus strikes a unique balance. On the one hand, we provide a comprehensive overview of the available modelling techniques and their applications, including key mathematical expressions. On the other hand, the content is accessible to mathematicians and biomedical scientists alike to accommodate the interdisciplinary nature of cancer research.

Influence of the nuclear membrane, active transport, and cell shape on the Hes1 and p53-Mdm2 pathways: insights from spatio-temporal modelling.

There are many intracellular signalling pathways where the spatial distribution of the molecular species cannot be neglected. These pathways often contain negative feedback loops and can exhibit oscillatory dynamics in space and time. Two such pathways are those involving Hes1 and p53-Mdm2, both of which are implicated in cancer. In this paper we further develop the partial differential equation (PDE) models of Sturrock et al. (J. Theor. Biol., 273:15-31, 2011) which were used to study these dynamics. We extend these PDE models by including a nuclear membrane and active transport, assuming that proteins are convected in the cytoplasm towards the nucleus in order to model transport along microtubules. We also account for Mdm2 inhibition of p53 transcriptional activity. Through numerical simulations we find ranges of values for the model parameters such that sustained oscillatory dynamics occur, consistent with available experimental measurements. We also find that our model extensions act to broaden the parameter ranges that yield oscillations. Hence oscillatory behaviour is made more robust by the inclusion of both the nuclear membrane and active transport. In order to bridge the gap between in vivo and in silico experiments, we investigate more realistic cell geometries by using an imported image of a real cell as our computational domain. For the extended p53-Mdm2 model, we consider the effect of microtubule-disrupting drugs and proteasome inhibitor drugs, obtaining results that are in agreement with experimental studies.

Fine-mapping and association analysis of candidate genes for papilla number in sea cucumber, Apostichopus japonicus.

The papilla number is one of the most economically important traits of sea cucumber in the China marketing trade. However, the genetic basis for papilla number diversity in holothurians is still scarce. In the present study, we conducted genome-wide association studies (GWAS) for the trait papilla number of sea cucumbers utilizing a set of 400,186 high-quality SNPs derived from 200 sea cucumbers. Two significant trait-associated SNPs that passed Bonferroni correction (P < 1.25E-7) were located in the intergenic region near PATS1 and the genic region of EIF4G, which were reported to play a pivotal role in cell growth and proliferation. The fine-mapping regions around the top two lead SNPs provided precise causative loci/genes related to papilla formation and cellular activity, including PPP2R3C, GBP1, and BCAS3. Potential SNPs with P < 1E-4 were acquired for the following GO and KEGG enrichment analysis. Moreover, the two lead SNPs were verified in another population of sea cucumber, and the expressive detection of three potential candidate genes PATS1, PPP2R3C, and EIF4G that near or cover the two lead SNPs was conducted in papilla tissue of TG (Top papilla number group) and BG (Bottom papilla number group) by qRT-PCR. We found the significantly higher expression profile of PATS1 (3.34-fold), PPP2R3C (4.90-fold), and EIF4G (4.23-fold) in TG, implying their potential function in papilla polymorphism. The present results provide valuable information to decipher the phenotype differences of the papilla trait and will provide a scientific basis for selective breeding in sea cucumbers. Supplementary Information: The online version contains supplementary material available at 10.1007/s42995-022-00139-w.

Emergent expression of fitness-conferring genes by phenotypic selection.

Genotypic and phenotypic adaptation is the consequence of ongoing natural selection in populations and is key to predicting and preventing drug resistance. Whereas classic antibiotic persistence is all-or-nothing, here we demonstrate that an antibiotic resistance gene displays linear dose-responsive selection for increased expression in proportion to rising antibiotic concentration in growing Escherichia coli populations. Furthermore, we report the potentially wide-spread nature of this form of emergent gene expression (EGE) by instantaneous phenotypic selection process under bactericidal and bacteriostatic antibiotic treatment, as well as an amino acid synthesis pathway enzyme under a range of auxotrophic conditions. We propose an analogy to Ohm's law in electricity (V = IR), where selection pressure acts similarly to voltage (V), gene expression to current (I), and resistance (R) to cellular machinery constraints and costs. Lastly, mathematical modeling using agent-based models of stochastic gene expression in growing populations and Bayesian model selection reveal that the EGE mechanism requires variability in gene expression within an isogenic population, and a cellular "memory" from positive feedbacks between growth and expression of any fitness-conferring gene. Finally, we discuss the connection of the observed phenomenon to a previously described general fluctuation-response relationship in biology.

Spatio-temporal modelling of the Hes1 and p53-Mdm2 intracellular signalling pathways.

The correct localisation of transcription factors is vitally important for the proper functioning of many intracellular signalling pathways. Experimental data has shown that many pathways exhibit oscillations in concentrations of the substances involved, both temporally and spatially. Negative feedback loops are important components of these oscillations, providing fine regulation for the factors involved. In this paper we consider mathematical models of two such pathways-Hes1 and p53-Mdm2. Building on previous mathematical modelling approaches, we derive systems of partial differential equations to capture the evolution in space and time of the variables in the Hes1 and p53-Mdm2 systems. Through computational simulations we show that our reaction-diffusion models are able to produce sustained oscillations both spatially and temporally, accurately reflecting experimental evidence and advancing previous models. The simulations of our models also allow us to calculate a diffusion coefficient range for the variables in each mRNA and protein system, as well as ranges for other key parameters of the models, where sustained oscillations are observed. Finally, by exploiting the explicitly spatial nature of the partial differential equations, we are also able to manipulate mathematically the spatial location of the ribosomes, thus controlling where the proteins are synthesized within the cytoplasm. The results of these simulations predict an optimal distance outside the nucleus where protein synthesis should take place in order to generate sustained oscillations. Using partial differential equation models, new information can be gained about the precise spatio-temporal dynamics of mRNA and proteins. The ability to determine spatial localisation of proteins within the cell is likely to yield fresh insight into a range of cellular diseases such as diabetes and cancer.

Phenotypic selection through cell death: stochastic modelling of O-6-methylguanine-DNA methyltransferase dynamics.

Glioblastoma (GBM) is the most aggressive malignant primary brain tumour with a median overall survival of 15 months. To treat GBM, patients currently undergo a surgical resection followed by exposure to radiotherapy and concurrent and adjuvant temozolomide (TMZ) chemotherapy. However, this protocol often leads to treatment failure, with drug resistance being the main reason behind this. To date, many studies highlight the role of O-6-methylguanine-DNA methyltransferase (MGMT) in conferring drug resistance. The mechanism through which MGMT confers resistance is not well studied-particularly in terms of computational models. With only a few reasonable biological assumptions, we were able to show that even a minimal model of MGMT expression could robustly explain TMZ-mediated drug resistance. In particular, we showed that for a wide range of parameter values constrained by novel cell growth and viability assays, a model accounting for only stochastic gene expression of MGMT coupled with cell growth, division, partitioning and death was able to exhibit phenotypic selection of GBM cells expressing MGMT in response to TMZ. Furthermore, we found this selection allowed the cells to pass their acquired phenotypic resistance onto daughter cells in a stable manner (as long as TMZ is provided). This suggests that stochastic gene expression alone is enough to explain the development of chemotherapeutic resistance.

Identifying ultrasensitive HGF dose-response functions in a 3D mammalian system for synthetic morphogenesis.

Nonlinear responses to signals are widespread natural phenomena that affect various cellular processes. Nonlinearity can be a desirable characteristic for engineering living organisms because it can lead to more switch-like responses, similar to those underlying the wiring in electronics. Steeper functions are described as ultrasensitive, and can be applied in synthetic biology by using various techniques including receptor decoys, multiple co-operative binding sites, and sequential positive feedbacks. Here, we explore the inherent non-linearity of a biological signaling system to identify functions that can potentially be exploited using cell genome engineering. For this, we performed genome-wide transcription profiling to identify genes with ultrasensitive response functions to Hepatocyte Growth Factor (HGF). We identified 3,527 genes that react to increasing concentrations of HGF, in Madin-Darby canine kidney (MDCK) cells, grown as cysts in 3D collagen cell culture. By fitting a generic Hill function to the dose-responses of these genes we obtained a measure of the ultrasensitivity of HGF-responsive genes, identifying a subset with higher apparent Hill coefficients (e.g. MMP1, TIMP1, SNORD75, SNORD86 and ERRFI1). The regulatory regions of these genes are potential candidates for future engineering of synthetic mammalian gene circuits requiring nonlinear responses to HGF signalling.

Protein abundance may regulate sensitivity to external cues in polarized cells.

Cell polarization is a ubiquitous process which results in cellular constituents being organized into discrete intracellular spatial domains. It occurs in a variety of cell types, including epithelial cells, immune system cells and neurons. A key player in this process is the Par protein family whose asymmetric localization to anterior and posterior parts of the cell is crucial for proper division and cell fate specification. In this paper, we explore a stochastic analogue of the temporal model of Par protein interactions first developed in Dawes & Munro (Dawes and Munro 2011 Biophys. J. 101, 1412-1422. (doi:10.1016/j.bpj.2011.07.030)). We focus on how protein abundance influences the behaviour of both the deterministic and stochastic versions of the model. In Dawes & Munro (2011), it was found that bistable behaviour in the temporal model of Par protein led to the existence of complementary domains in the corresponding spatio-temporal model. Here, we find that the corresponding temporal stochastic model permits switching behaviour (the model solution 'jumps' between steady states) for lower protein abundances, whereas for higher protein abundances the stochastic and deterministic models are in good agreement (the model solution evolves to one of two steady states). This led us to the testable hypothesis that cells with lower abundances of Par protein may be more sensitive to external cues, whereas cells with higher abundances of Par protein may be less sensitive to external cues. In order to gain more control over the precise abundance of Par protein, we proposed and explored a second model (again, examining both deterministic and stochastic versions) in which the total number of Par molecules is conserved. We found that this model required an additional dimerization reaction in the cytoplasm in order for bistable and switching behaviour to be found. Once this additional reaction was included, we found that both the first and second models gave qualitatively similar results but in different regions of the parameter space, suggesting a further regulatory mechanism that cells could potentially use to modulate their response to external signals.

Spatial stochastic modelling of the Hes1 gene regulatory network: intrinsic noise can explain heterogeneity in embryonic stem cell differentiation.

Individual mouse embryonic stem cells have been found to exhibit highly variable differentiation responses under the same environmental conditions. The noisy cyclic expression of Hes1 and its downstream genes are known to be responsible for this, but the mechanism underlying this variability in expression is not well understood. In this paper, we show that the observed experimental data and diverse differentiation responses can be explained by a spatial stochastic model of the Hes1 gene regulatory network. We also propose experiments to control the precise differentiation response using drug treatment.

A spatio-temporal model of Notch signalling in the zebrafish segmentation clock: conditions for synchronised oscillatory dynamics.

In the vertebrate embryo, tissue blocks called somites are laid down in head-to-tail succession, a process known as somitogenesis. Research into somitogenesis has been both experimental and mathematical. For zebrafish, there is experimental evidence for oscillatory gene expression in cells in the presomitic mesoderm (PSM) as well as evidence that Notch signalling synchronises the oscillations in neighbouring PSM cells. A biological mechanism has previously been proposed to explain these phenomena. Here we have converted this mechanism into a mathematical model of partial differential equations in which the nuclear and cytoplasmic diffusion of protein and mRNA molecules is explicitly considered. By performing simulations, we have found ranges of values for the model parameters (such as diffusion and degradation rates) that yield oscillatory dynamics within PSM cells and that enable Notch signalling to synchronise the oscillations in two touching cells. Our model contains a Hill coefficient that measures the co-operativity between two proteins (Her1, Her7) and three genes (her1, her7, deltaC) which they inhibit. This coefficient appears to be bounded below by the requirement for oscillations in individual cells and bounded above by the requirement for synchronisation. Consistent with experimental data and a previous spatially non-explicit mathematical model, we have found that signalling can increase the average level of Her1 protein. Biological pattern formation would be impossible without a certain robustness to variety in cell shape and size; our results possess such robustness. Our spatially-explicit modelling approach, together with new imaging technologies that can measure intracellular protein diffusion rates, is likely to yield significant new insight into somitogenesis and other biological processes.