Parameter estimation and inference for stochastic reaction-diffusion systems: application to morphogenesis in D. melanogaster

Inferring Stochastic Rates from Heterogeneous Snapshots of Particle Positions

Many imaging techniques for biological systems-like fixation of cells coupled with fluorescence microscopy-provide sharp spatial resolution in reporting locations of individuals at a single moment in time but also destroy the dynamics they intend to capture. These snapshot observations contain no information about individual trajectories, but still encode information about movement and demographic dynamics, especially when combined with a well-motivated biophysical model. The relationship between spatially evolving populations and single-moment representations of their collective locations is well-established with partial differential equations (PDEs) and their inverse problems. However, experimental data is commonly a set of locations whose number is insufficient to approximate a continuous-in-space PDE solution. Here, motivated by popular subcellular imaging data of gene expression, we embrace the stochastic nature of the data and investigate the mathematical foundations of parametrically inferring demographic rates from snapshots of particles undergoing birth, diffusion, and death in a nuclear or cellular domain. Toward inference, we rigorously derive a connection between individual particle paths and their presentation as a Poisson spatial process. Using this framework, we investigate the properties of the resulting inverse problem and study factors that affect quality of inference. One pervasive feature of this experimental regime is the presence of cell-to-cell heterogeneity. Rather than being a hindrance, we show that cell-to-cell geometric heterogeneity can increase the quality of inference on dynamics for certain parameter regimes. Altogether, the results serve as a basis for more detailed investigations of subcellular spatial patterns of RNA molecules and other stochastically evolving populations that can only be observed for single instants in their time evolution.

Inferring stochastic rates from heterogeneous snapshots of particle positions

Many imaging techniques for biological systems - like fixation of cells coupled with fluorescence microscopy - provide sharp spatial resolution in reporting locations of individuals at a single moment in time but also destroy the dynamics they intend to capture. These snapshot observations contain no information about individual trajectories, but still encode information about movement and demographic dynamics, especially when combined with a well-motivated biophysical model. The relationship between spatially evolving populations and single-moment representations of their collective locations is well-established with partial differential equations (PDEs) and their inverse problems. However, experimental data is commonly a set of locations whose number is insufficient to approximate a continuous-in-space PDE solution. Here, motivated by popular subcellular imaging data of gene expression, we embrace the stochastic nature of the data and investigate the mathematical foundations of parametrically inferring demographic rates from snapshots of particles undergoing birth, diffusion, and death in a nuclear or cellular domain. Toward inference, we rigorously derive a connection between individual particle paths and their presentation as a Poisson spatial process. Using this framework, we investigate the properties of the resulting inverse problem and study factors that affect quality of inference. One pervasive feature of this experimental regime is the presence of cell-to-cell heterogeneity. Rather than being a hindrance, we show that cell-to-cell geometric heterogeneity can increase the quality of inference on dynamics for certain parameter regimes. Altogether, the results serve as a basis for more detailed investigations of subcellular spatial patterns of RNA molecules and other stochastically evolving populations that can only be observed for single instants in their time evolution.

Bespoke Turing Systems

Reaction-diffusion systems are an intensively studied form of partial differential equation, frequently used to produce spatially heterogeneous patterned states from homogeneous symmetry breaking via the Turing instability. Although there are many prototypical "Turing systems" available, determining their parameters, functional forms, and general appropriateness for a given application is often difficult. Here, we consider the reverse problem. Namely, suppose we know the parameter region associated with the reaction kinetics in which patterning is required-we present a constructive framework for identifying systems that will exhibit the Turing instability within this region, whilst in addition often allowing selection of desired patterning features, such as spots, or stripes. In particular, we show how to build a system of two populations governed by polynomial morphogen kinetics such that the: patterning parameter domain (in any spatial dimension), morphogen phases (in any spatial dimension), and even type of resulting pattern (in up to two spatial dimensions) can all be determined. Finally, by employing spatial and temporal heterogeneity, we demonstrate that mixed mode patterns (spots, stripes, and complex prepatterns) are also possible, allowing one to build arbitrarily complicated patterning landscapes. Such a framework can be employed pedagogically, or in a variety of contemporary applications in designing synthetic chemical and biological patterning systems. We also discuss the implications that this freedom of design has on using reaction-diffusion systems in biological modelling and suggest that stronger constraints are needed when linking theory and experiment, as many simple patterns can be easily generated given freedom to choose reaction kinetics.

Cox process representation and inference for stochastic reaction-diffusion processes

Complex behaviour in many systems arises from the stochastic interactions of spatially distributed particles or agents. Stochastic reaction–diffusion processes are widely used to model such behaviour in disciplines ranging from biology to the social sciences, yet they are notoriously difficult to simulate and calibrate to observational data. Here we use ideas from statistical physics and machine learning to provide a solution to the inverse problem of learning a stochastic reaction–diffusion process from data. Our solution relies on a non-trivial connection between stochastic reaction–diffusion processes and spatio-temporal Cox processes, a well-studied class of models from computational statistics. This connection leads to an efficient and flexible algorithm for parameter inference and model selection. Our approach shows excellent accuracy on numeric and real data examples from systems biology and epidemiology. Our work provides both insights into spatio-temporal stochastic systems, and a practical solution to a long-standing problem in computational modelling.

Stochastic reaction-diffusion systems are used for modelling spatial dynamics in many disciplines, but parameter inference and model selection remain challenging. Here the authors offer a solution enabled by a connection between reaction-diffusion and the well-studied spatio-temporal Cox processes.

Parameter identification problems in the modelling of cell motility

We present a novel parameter identification algorithm for the estimation of parameters in models of cell motility using imaging data of migrating cells. Two alternative formulations of the objective functional that measures the difference between the computed and observed data are proposed and the parameter identification problem is formulated as a minimisation problem of nonlinear least squares type. A Levenberg-Marquardt based optimisation method is applied to the solution of the minimisation problem and the details of the implementation are discussed. A number of numerical experiments are presented which illustrate the robustness of the algorithm to parameter identification in the presence of large deformations and noisy data and parameter identification in three dimensional models of cell motility. An application to experimental data is also presented in which we seek to identify parameters in a model for the monopolar growth of fission yeast cells using experimental imaging data. Our numerical tests allow us to compare the method with the two different formulations of the objective functional and we conclude that the results with both objective functionals seem to agree.

Bayesian inference of the lung alveolar spatial model for the identification of alveolar mechanics associated with acute respiratory distress syndrome

Acute respiratory distress syndrome (ARDS) is acute lung failure secondary to severe systemic inflammation, resulting in a derangement of alveolar mechanics (i.e. the dynamic change in alveolar size and shape during tidal ventilation), leading to alveolar instability that can cause further damage to the pulmonary parenchyma. Mechanical ventilation is a mainstay in the treatment of ARDS, but may induce mechano-physical stresses on unstable alveoli, which can paradoxically propagate the cellular and molecular processes exacerbating ARDS pathology. This phenomenon is called ventilator induced lung injury (VILI), and plays a significant role in morbidity and mortality associated with ARDS. In order to identify optimal ventilation strategies to limit VILI and treat ARDS, it is necessary to understand the complex interplay between biological and physical mechanisms of VILI, first at the alveolar level, and then in aggregate at the whole-lung level. Since there is no current consensus about the underlying dynamics of alveolar mechanics, as an initial step we investigate the ventilatory dynamics of an alveolar sac (AS) with the lung alveolar spatial model (LASM), a 3D spatial biomechanical representation of the AS and its interaction with airflow pressure and the surface tension effects of pulmonary surfactant. We use the LASM to identify the mechanical ramifications of alveolar dynamics associated with ARDS. Using graphical processing unit parallel algorithms, we perform Bayesian inference on the model parameters using experimental data from rat lung under control and Tween-induced ARDS conditions. Our results provide two plausible models that recapitulate two fundamental hypotheses about volume change at the alveolar level: (1) increase in alveolar size through isotropic volume change, or (2) minimal change in AS radius with primary expansion of the mouth of the AS, with the implication that the majority of change in lung volume during the respiratory cycle occurs in the alveolar ducts. These two model solutions correspond to significantly different mechanical properties of the tissue, and we discuss the implications of these different properties and the requirements for new experimental data to discriminate between the hypotheses.

Efficient reverse-engineering of a developmental gene regulatory network

Understanding the complex regulatory networks underlying development and evolution of multi-cellular organisms is a major problem in biology. Computational models can be used as tools to extract the regulatory structure and dynamics of such networks from gene expression data. This approach is called reverse engineering. It has been successfully applied to many gene networks in various biological systems. However, to reconstitute the structure and non-linear dynamics of a developmental gene network in its spatial context remains a considerable challenge. Here, we address this challenge using a case study: the gap gene network involved in segment determination during early development of Drosophila melanogaster. A major problem for reverse-engineering pattern-forming networks is the significant amount of time and effort required to acquire and quantify spatial gene expression data. We have developed a simplified data processing pipeline that considerably increases the throughput of the method, but results in data of reduced accuracy compared to those previously used for gap gene network inference. We demonstrate that we can infer the correct network structure using our reduced data set, and investigate minimal data requirements for successful reverse engineering. Our results show that timing and position of expression domain boundaries are the crucial features for determining regulatory network structure from data, while it is less important to precisely measure expression levels. Based on this, we define minimal data requirements for gap gene network inference. Our results demonstrate the feasibility of reverse-engineering with much reduced experimental effort. This enables more widespread use of the method in different developmental contexts and organisms. Such systematic application of data-driven models to real-world networks has enormous potential. Only the quantitative investigation of a large number of developmental gene regulatory networks will allow us to discover whether there are rules or regularities governing development and evolution of complex multi-cellular organisms.

State and parameter estimation of the heat shock response system using Kalman and particle filters

MOTIVATION

Traditional models of systems biology describe dynamic biological phenomena as solutions to ordinary differential equations, which, when parameters in them are set to correct values, faithfully mimic observations. Often parameter values are tweaked by hand until desired results are achieved, or computed from biochemical experiments carried out in vitro. Of interest in this article, is the use of probabilistic modelling tools with which parameters and unobserved variables, modelled as hidden states, can be estimated from limited noisy observations of parts of a dynamical system.

RESULTS

Here we focus on sequential filtering methods and take a detailed look at the capabilities of three members of this family: (i) extended Kalman filter (EKF), (ii) unscented Kalman filter (UKF) and (iii) the particle filter, in estimating parameters and unobserved states of cellular response to sudden temperature elevation of the bacterium Escherichia coli. While previous literature has studied this system with the EKF, we show that parameter estimation is only possible with this method when the initial guesses are sufficiently close to the true values. The same turns out to be true for the UKF. In this thorough empirical exploration, we show that the non-parametric method of particle filtering is able to reliably estimate parameters and states, converging from initial distributions relatively far away from the underlying true values.

AVAILABILITY AND IMPLEMENTATION

Software implementation of the three filters on this problem can be freely downloaded from http://users.ecs.soton.ac.uk/mn/HeatShock

The neutrophil's eye-view: inference and visualisation of the chemoattractant field driving cell chemotaxis in vivo

As we begin to understand the signals that drive chemotaxis in vivo, it is becoming clear that there is a complex interplay of chemotactic factors, which changes over time as the inflammatory response evolves. New animal models such as transgenic lines of zebrafish, which are near transparent and where the neutrophils express a green fluorescent protein, have the potential to greatly increase our understanding of the chemotactic process under conditions of wounding and infection from video microscopy data. Measurement of the chemoattractants over space (and their evolution over time) is a key objective for understanding the signals driving neutrophil chemotaxis. However, it is not possible to measure and visualise the most important contributors to in vivo chemotaxis, and in fact the understanding of the main contributors at any particular time is incomplete. The key insight that we make in this investigation is that the neutrophils themselves are sensing the underlying field that is driving their action and we can use the observations of neutrophil movement to infer the hidden net chemoattractant field by use of a novel computational framework. We apply the methodology to multiple in vivo neutrophil recruitment data sets to demonstrate this new technique and find that the method provides consistent estimates of the chemoattractant field across the majority of experiments. The framework that we derive represents an important new methodology for cell biologists investigating the signalling processes driving cell chemotaxis, which we label the neutrophils eye-view of the chemoattractant field.

Gaussian process modelling for bicoid mRNA regulation in spatio-temporal Bicoid profile

MOTIVATION

Bicoid protein molecules, translated from maternally provided bicoid mRNA, establish a concentration gradient in Drosophila early embryonic development. There is experimental evidence that the synthesis and subsequent destruction of this protein is regulated at source by precise control of the stability of the maternal mRNA. Can we infer the driving function at the source from noisy observations of the spatio-temporal protein profile? We use non-parametric Gaussian process regression for modelling the propagation of Bicoid in the embryo and infer aspects of source regulation as a posterior function.

RESULTS

With synthetic data from a 1D diffusion model with a source simulated to model mRNA stability regulation, our results establish that the Gaussian process method can accurately infer the driving function and capture the spatio-temporal dynamics of embryonic Bicoid propagation. On real data from the FlyEx database, too, the reconstructed source function is indicative of stability regulation, but is temporally smoother than what we expected, partly due to the fact that the dataset is only partially observed. To be in line with recent thinking on the subject, we also analyse this model with a spatial gradient of maternal mRNA, rather than being fixed at only the anterior pole.

CONTACT

m.niranjan@southampton.ac.uk

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

The role of regulated mRNA stability in establishing bicoid morphogen gradient in Drosophila embryonic development

The Bicoid morphogen is amongst the earliest triggers of differential spatial pattern of gene expression and subsequent cell fate determination in the embryonic development of Drosophila. This maternally deposited morphogen is thought to diffuse in the embryo, establishing a concentration gradient which is sensed by downstream genes. In most model based analyses of this process, the translation of the bicoid mRNA is thought to take place at a fixed rate from the anterior pole of the embryo and a supply of the resulting protein at a constant rate is assumed. Is this process of morphogen generation a passive one as assumed in the modelling literature so far, or would available data support an alternate hypothesis that the stability of the mRNA is regulated by active processes? We introduce a model in which the stability of the maternal mRNA is regulated by being held constant for a length of time, followed by rapid degradation. With this more realistic model of the source, we have analysed three computational models of spatial morphogen propagation along the anterior-posterior axis: (a) passive diffusion modelled as a deterministic differential equation, (b) diffusion enhanced by a cytoplasmic flow term; and (c) diffusion modelled by stochastic simulation of the corresponding chemical reactions. Parameter estimation on these models by matching to publicly available data on spatio-temporal Bicoid profiles suggests strong support for regulated stability over either a constant supply rate or one where the maternal mRNA is permitted to degrade in a passive manner.