Contributions of cell behavior to geometric order in embryonic cartilage.

During early development, cartilage provides shape and stability to the embryo while serving as a precursor for the skeleton. Correct formation of embryonic cartilage is hence essential for healthy development. In vertebrate cranial cartilage, it has been observed that a flat and laterally extended macroscopic geometry is linked to regular microscopic structure consisting of tightly packed, short, transversal clonar columns. However, it remains an ongoing challenge to identify how individual cells coordinate to successfully shape the tissue, and more precisely which mechanical interactions and cell behaviors contribute to the generation and maintenance of this columnar cartilage geometry during embryogenesis. Here, we apply a three-dimensional cell-based computational model to investigate mechanical principles contributing to column formation. The model accounts for clonal expansion, anisotropic proliferation and the geometrical arrangement of progenitor cells in space. We confirm that oriented cell divisions and repulsive mechanical interactions between cells are key drivers of column formation. In addition, the model suggests that column formation benefits from the spatial gaps created by the extracellular matrix in the initial configuration, and that column maintenance is facilitated by sequential proliferative phases. Our model thus correctly predicts the dependence of local order on division orientation and tissue thickness. The present study presents the first cell-based simulations of cell mechanics during cranial cartilage formation and we anticipate that it will be useful in future studies on the formation and growth of other cartilage geometries.

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.

Identification of dynamic mass-action biochemical reaction networks using sparse Bayesian methods.

Identifying the reactions that govern a dynamical biological system is a crucial but challenging task in systems biology. In this work, we present a data-driven method to infer the underlying biochemical reaction system governing a set of observed species concentrations over time. We formulate the problem as a regression over a large, but limited, mass-action constrained reaction space and utilize sparse Bayesian inference via the regularized horseshoe prior to produce robust, interpretable biochemical reaction networks, along with uncertainty estimates of parameters. The resulting systems of chemical reactions and posteriors inform the biologist of potentially several reaction systems that can be further investigated. We demonstrate the method on two examples of recovering the dynamics of an unknown reaction system, to illustrate the benefits of improved accuracy and information obtained.

CBMOS: a GPU-enabled Python framework for the numerical study of center-based models.

BACKGROUND: Cell-based models are becoming increasingly popular for applications in developmental biology. However, the impact of numerical choices on the accuracy and efficiency of the simulation of these models is rarely meticulously tested. Without concrete studies to differentiate between solid model conclusions and numerical artifacts, modelers are at risk of being misled by their experiments' results. Most cell-based modeling frameworks offer a feature-rich environment, providing a wide range of biological components, but are less suitable for numerical studies. There is thus a need for software specifically targeted at this use case. RESULTS: We present CBMOS, a Python framework for the simulation of the center-based or cell-centered model. Contrary to other implementations, CBMOS' focus is on facilitating numerical study of center-based models by providing access to multiple ordinary differential equation solvers and force functions through a flexible, user-friendly interface and by enabling rapid testing through graphics processing unit (GPU) acceleration. We show-case its potential by illustrating two common workflows: (1) comparison of the numerical properties of two solvers within a Jupyter notebook and (2) measuring average wall times of both solvers on a high performance computing cluster. More specifically, we confirm that although for moderate accuracy levels the backward Euler method allows for larger time step sizes than the commonly used forward Euler method, its additional computational cost due to being an implicit method prohibits its use for practical test cases. CONCLUSIONS: CBMOS is a flexible, easy-to-use Python implementation of the center-based model, exposing both basic model assumptions and numerical components to the user. It is available on GitHub and PyPI under an MIT license. CBMOS allows for fast prototyping on a central processing unit for small systems through the use of NumPy. Using CuPy on a GPU, cell populations of up to 10,000 cells can be simulated within a few seconds. As such, it will substantially lower the time investment for any modeler to check the crucial assumption that model conclusions are independent of numerical issues.

Scalable machine learning-assisted model exploration and inference using Sciope.

SUMMARY: Discrete stochastic models of gene regulatory networks are fundamental tools for in silico study of stochastic gene regulatory networks. Likelihood-free inference and model exploration are critical applications to study a system using such models. However, the massive computational cost of complex, high-dimensional and stochastic modelling currently limits systematic investigation to relatively simple systems. Recently, machine-learning-assisted methods have shown great promise to handle larger, more complex models. To support both ease-of-use of this new class of methods, as well as their further development, we have developed the scalable inference, optimization and parameter exploration (Sciope) toolbox. Sciope is designed to support new algorithms for machine-learning-assisted model exploration and likelihood-free inference. Moreover, it is built ground up to easily leverage distributed and heterogeneous computational resources for convenient parallelism across platforms from workstations to clouds. AVAILABILITY AND IMPLEMENTATION: The Sciope Python3 toolbox is freely available on https://github.com/Sciope/Sciope, and has been tested on Linux, Windows and macOS platforms. SUPPLEMENTARY INFORMATION: Supplementary information is available at Bioinformatics online.

Epidemiological modeling in StochSS Live!

SUMMARY: We present StochSS Live!, a web-based service for modeling, simulation and analysis of a wide range of mathematical, biological and biochemical systems. Using an epidemiological model of COVID-19, we demonstrate the power of StochSS Live! to enable researchers to quickly develop a deterministic or a discrete stochastic model, infer its parameters and analyze the results. AVAILABILITY AND IMPLEMENTATION: StochSS Live! is freely available at https://live.stochss.org/. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Accelerated regression-based summary statistics for discrete stochastic systems via approximate simulators.

BACKGROUND: Approximate Bayesian Computation (ABC) has become a key tool for calibrating the parameters of discrete stochastic biochemical models. For higher dimensional models and data, its performance is strongly dependent on having a representative set of summary statistics. While regression-based methods have been demonstrated to allow for the automatic construction of effective summary statistics, their reliance on first simulating a large training set creates a significant overhead when applying these methods to discrete stochastic models for which simulation is relatively expensive. In this τ work, we present a method to reduce this computational burden by leveraging approximate simulators of these systems, such as ordinary differential equations and τ-Leaping approximations. RESULTS: We have developed an algorithm to accelerate the construction of regression-based summary statistics for Approximate Bayesian Computation by selectively using the faster approximate algorithms for simulations. By posing the problem as one of ratio estimation, we use state-of-the-art methods in machine learning to show that, in many cases, our algorithm can significantly reduce the number of simulations from the full resolution model at a minimal cost to accuracy and little additional tuning from the user. We demonstrate the usefulness and robustness of our method with four different experiments. CONCLUSIONS: We provide a novel algorithm for accelerating the construction of summary statistics for stochastic biochemical systems. Compared to the standard practice of exclusively training from exact simulator samples, our method is able to dramatically reduce the number of required calls to the stochastic simulator at a minimal loss in accuracy. This can immediately be implemented to increase the overall speed of the ABC workflow for estimating parameters in complex systems.

A multiscale compartment-based model of stochastic gene regulatory networks using hitting-time analysis.

Spatial stochastic models of single cell kinetics are capable of capturing both fluctuations in molecular numbers and the spatial dependencies of the key steps of intracellular regulatory networks. The spatial stochastic model can be simulated both on a detailed microscopic level using particle tracking and on a mesoscopic level using the reaction-diffusion master equation. However, despite substantial progress on simulation efficiency for spatial models in the last years, the computational cost quickly becomes prohibitively expensive for tasks that require repeated simulation of thousands or millions of realizations of the model. This limits the use of spatial models in applications such as multicellular simulations, likelihood-free parameter inference, and robustness analysis. Further approximation of the spatial dynamics is needed to accelerate such computational engineering tasks. We here propose a multiscale model where a compartment-based model approximates a detailed spatial stochastic model. The compartment model is constructed via a first-exit time analysis on the spatial model, thus capturing critical spatial aspects of the fine-grained simulations, at a cost close to the simple well-mixed model. We apply the multiscale model to a canonical model of negative-feedback gene regulation, assess its accuracy over a range of parameters, and demonstrate that the approximation can yield substantial speedups for likelihood-free parameter inference.

Smart computational exploration of stochastic gene regulatory network models using human-in-the-loop semi-supervised learning.

MOTIVATION: Discrete stochastic models of gene regulatory network models are indispensable tools for biological inquiry since they allow the modeler to predict how molecular interactions give rise to nonlinear system output. Model exploration with the objective of generating qualitative hypotheses about the workings of a pathway is usually the first step in the modeling process. It involves simulating the gene network model under a very large range of conditions, due to the large uncertainty in interactions and kinetic parameters. This makes model exploration highly computational demanding. Furthermore, with no prior information about the model behavior, labor-intensive manual inspection of very large amounts of simulation results becomes necessary. This limits systematic computational exploration to simplistic models. RESULTS: We have developed an interactive, smart workflow for model exploration based on semi-supervised learning and human-in-the-loop labeling of data. The workflow lets a modeler rapidly discover ranges of interesting behaviors predicted by the model. Utilizing that similar simulation output is in proximity of each other in a feature space, the modeler can focus on informing the system about what behaviors are more interesting than others by labeling, rather than analyzing simulation results with custom scripts and workflows. This results in a large reduction in time-consuming manual work by the modeler early in a modeling project, which can substantially reduce the time needed to go from an initial model to testable predictions and downstream analysis. AVAILABILITY AND IMPLEMENTATION: A python-package is available at https://github.com/Wrede/mio.git. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Impact of Force Function Formulations on the Numerical Simulation of Centre-Based Models.

Centre-based or cell-centre models are a framework for the computational study of multicellular systems with widespread use in cancer modelling and computational developmental biology. At the core of these models are the numerical method used to update cell positions and the force functions that encode the pairwise mechanical interactions of cells. For the latter, there are multiple choices that could potentially affect both the biological behaviour captured, and the robustness and efficiency of simulation. For example, available open-source software implementations of centre-based models rely on different force functions for their default behaviour and it is not straightforward for a modeller to know if these are interchangeable. Our study addresses this problem and contributes to the understanding of the potential and limitations of three popular force functions from a numerical perspective. We show empirically that choosing the force parameters such that the relaxation time for two cells after cell division is consistent between different force functions results in good agreement of the population radius of a two-dimensional monolayer relaxing mechanically after intense cell proliferation. Furthermore, we report that numerical stability is not sufficient to prevent unphysical cell trajectories following cell division, and consequently, that too large time steps can cause geometrical differences at the population level.

Hierarchical algorithm for the reaction-diffusion master equation.

We have developed an algorithm coupling mesoscopic simulations on different levels in a hierarchy of Cartesian meshes. Based on the multiscale nature of the chemical reactions, some molecules in the system will live on a fine-grained mesh, while others live on a coarse-grained mesh. By allowing molecules to transfer from the fine levels to the coarse levels when appropriate, we show that we can save up to three orders of magnitude of computational time compared to microscopic simulations or highly resolved mesoscopic simulations, without losing significant accuracy. We demonstrate this in several numerical examples with systems that cannot be accurately simulated with a coarse-grained mesoscopic model.

A 3D Multiscale Model to Explore the Role of EGFR Overexpression in Tumourigenesis.

The epidermal growth factor receptor (EGFR) signalling cascade is one of the main pathways that regulate the survival and division of mammalian cells. It is also one of the most altered transduction pathways in cancer. Acquired mutations in the EGFR/ERK pathway can cause the overexpression of EGFR on the surface of the cell, while others downregulate the inactivation of switched on intracellular proteins such as Ras and Raf. This upregulates the activity of ERK and promotes cell division. We develop a 3D multiscale model to explore the role of EGFR overexpression on tumour initiation. In this model, cells are described as individual objects that move, interact, divide, proliferate, and die by apoptosis. We use Brownian Dynamics to describe the extracellular and intracellular regulations of cells as well as the spatial and stochastic effects influencing them. The fate of each cell depends on the number of active transcription factors in the nucleus. We use numerical simulations to investigate the individual and combined effects of mutations on the intracellular regulation of individual cells. Next, we show that the distance between active receptors increase the level of EGFR/ERK signalling. We demonstrate the usefulness of the model by quantifying the impact of mutational alterations in the EGFR/ERK pathway on the growth rate of in silico tumours.

BAMSI: a multi-cloud service for scalable distributed filtering of massive genome data.

BACKGROUND: The advent of next-generation sequencing (NGS) has made whole-genome sequencing of cohorts of individuals a reality. Primary datasets of raw or aligned reads of this sort can get very large. For scientific questions where curated called variants are not sufficient, the sheer size of the datasets makes analysis prohibitively expensive. In order to make re-analysis of such data feasible without the need to have access to a large-scale computing facility, we have developed a highly scalable, storage-agnostic framework, an associated API and an easy-to-use web user interface to execute custom filters on large genomic datasets. RESULTS: We present BAMSI, a Software as-a Service (SaaS) solution for filtering of the 1000 Genomes phase 3 set of aligned reads, with the possibility of extension and customization to other sets of files. Unique to our solution is the capability of simultaneously utilizing many different mirrors of the data to increase the speed of the analysis. In particular, if the data is available in private or public clouds - an increasingly common scenario for both academic and commercial cloud providers - our framework allows for seamless deployment of filtering workers close to data. We show results indicating that such a setup improves the horizontal scalability of the system, and present a possible use case of the framework by performing an analysis of structural variation in the 1000 Genomes data set. CONCLUSIONS: BAMSI constitutes a framework for efficient filtering of large genomic data sets that is flexible in the use of compute as well as storage resources. The data resulting from the filter is assumed to be greatly reduced in size, and can easily be downloaded or routed into e.g. a Hadoop cluster for subsequent interactive analysis using Hive, Spark or similar tools. In this respect, our framework also suggests a general model for making very large datasets of high scientific value more accessible by offering the possibility for organizations to share the cost of hosting data on hot storage, without compromising the scalability of downstream analysis.

Stochastic Simulation Service: Bridging the Gap between the Computational Expert and the Biologist.

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.

Mesoscopic-microscopic spatial stochastic simulation with automatic system partitioning.

The reaction-diffusion master equation (RDME) is a model that allows for efficient on-lattice simulation of spatially resolved stochastic chemical kinetics. Compared to off-lattice hard-sphere simulations with Brownian dynamics or Green's function reaction dynamics, the RDME can be orders of magnitude faster if the lattice spacing can be chosen coarse enough. However, strongly diffusion-controlled reactions mandate a very fine mesh resolution for acceptable accuracy. It is common that reactions in the same model differ in their degree of diffusion control and therefore require different degrees of mesh resolution. This renders mesoscopic simulation inefficient for systems with multiscale properties. Mesoscopic-microscopic hybrid methods address this problem by resolving the most challenging reactions with a microscale, off-lattice simulation. However, all methods to date require manual partitioning of a system, effectively limiting their usefulness as "black-box" simulation codes. In this paper, we propose a hybrid simulation algorithm with automatic system partitioning based on indirect a priori error estimates. We demonstrate the accuracy and efficiency of the method on models of diffusion-controlled networks in 3D.

MESOSCOPIC MODELING OF STOCHASTIC REACTION-DIFFUSION KINETICS IN THE SUBDIFFUSIVE REGIME.

Subdiffusion has been proposed as an explanation of various kinetic phenomena inside living cells. In order to fascilitate large-scale computational studies of subdiffusive chemical processes, we extend a recently suggested mesoscopic model of subdiffusion into an accurate and consistent reaction-subdiffusion computational framework. Two different possible models of chemical reaction are revealed and some basic dynamic properties are derived. In certain cases those mesoscopic models have a direct interpretation at the macroscopic level as fractional partial differential equations in a bounded time interval. Through analysis and numerical experiments we estimate the macroscopic effects of reactions under subdiffusive mixing. The models display properties observed also in experiments: for a short time interval the behavior of the diffusion and the reaction is ordinary, in an intermediate interval the behavior is anomalous, and at long times the behavior is ordinary again.

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.

Perspective: Stochastic algorithms for chemical kinetics.

We outline our perspective on stochastic chemical kinetics, paying particular attention to numerical simulation algorithms. We first focus on dilute, well-mixed systems, whose description using ordinary differential equations has served as the basis for traditional chemical kinetics for the past 150 years. For such systems, we review the physical and mathematical rationale for a discrete-stochastic approach, and for the approximations that need to be made in order to regain the traditional continuous-deterministic description. We next take note of some of the more promising strategies for dealing stochastically with stiff systems, rare events, and sensitivity analysis. Finally, we review some recent efforts to adapt and extend the discrete-stochastic approach to systems that are not well-mixed. In that currently developing area, we focus mainly on the strategy of subdividing the system into well-mixed subvolumes, and then simulating diffusional transfers of reactant molecules between adjacent subvolumes together with chemical reactions inside the subvolumes.

Convolutional Neural Networks as Summary Statistics for Approximate Bayesian Computation.

Approximate Bayesian Computation is widely used in systems biology for inferring parameters in stochastic gene regulatory network models. Its performance hinges critically on the ability to summarize high-dimensional system responses such as time series into a few informative, low-dimensional summary statistics. The quality of those statistics acutely impacts the accuracy of the inference task. Existing methods to select the best subset out of a pool of candidate statistics do not scale well with large pools of several tens to hundreds of candidate statistics. Since high quality statistics are imperative for good performance, this becomes a serious bottleneck when performing inference on complex and high-dimensional problems. This paper proposes a convolutional neural network architecture for automatically learning informative summary statistics of temporal responses. We show that the proposed network can effectively circumvent the statistics selection problem of the preprocessing step for ABC inference. The proposed approach is demonstrated on two benchmark problem and one challenging inference problem learning parameters in a high-dimensional stochastic genetic oscillator. We also study the impact of experimental design on network performance by comparing different data richness and data acquisition strategies.

Living in darkness: Exploring adaptation of Proteus anguinus in 3 dimensions by X-ray imaging.

BACKGROUND: Lightless caves can harbour a wide range of living organisms. Cave animals have evolved a set of morphological, physiological, and behavioural adaptations known as troglomorphisms, enabling their survival in the perpetual darkness, narrow temperature and humidity ranges, and nutrient scarcity of the subterranean environment. In this study, we focused on adaptations of skull shape and sensory systems in the blind cave salamander, Proteus anguinus, also known as olm or simply proteus-the largest cave tetrapod and the only European amphibian living exclusively in subterranean environments. This extraordinary amphibian compensates for the loss of sight by enhanced non-visual sensory systems including mechanoreceptors, electroreceptors, and chemoreceptors. We compared developmental stages of P. anguinus with Ambystoma mexicanum, also known as axolotl, to make an exemplary comparison between cave- and surface-dwelling paedomorphic salamanders. FINDINGS: We used contrast-enhanced X-ray computed microtomography for the 3D segmentation of the soft tissues in the head of P. anguinus and A. mexicanum. Sensory organs were visualized to elucidate how the animal is adapted to living in complete darkness. X-ray microCT datasets were provided along with 3D models for larval, juvenile, and adult specimens, showing the cartilage of the chondrocranium and the position, shape, and size of the brain, eyes, and olfactory epithelium. CONCLUSIONS: P. anguinus still keeps some of its secrets. Our high-resolution X-ray microCT scans together with 3D models of the anatomical structures in the head may help to elucidate the nature and origin of the mechanisms behind its adaptations to the subterranean environment, which led to a series of troglomorphisms.