Parameter estimation and change-point detection from Dynamic Contrast Enhanced MRI data using stochastic differential equations

Dynamic Contrast Enhanced imaging (DCE-imaging) following a contrast agent bolus allows the extraction of information on tissue micro-vascularization. The dynamic signals obtained from DCE-imaging are modeled by pharmacokinetic compartmental models which integrate the Arterial Input Function. These models use ordinary differential equations (ODEs) to describe the exchanges between the arterial and capillary plasma and the extravascular-extracellular space. Their least squares fitting takes into account measurement noises but fails to deal with unpredictable fluctuations due to external/internal sources of variations (patients’ anxiety, time-varying parameters, measurement errors in the input function, etc.). Adding Brownian components to the ODEs leads to stochastic differential equations (SDEs). In DCE-imaging, SDEs are discretely observed with an additional measurement noise. We propose to estimate the parameters of these noisy SDEs by maximum likelihood, using the Kalman filter. In DCE-imaging, the contrast agent injected in vein arrives in plasma with an unknown time delay. The delay parameter induces a change-point in the drift of the SDE and ODE models, which is estimated also. Estimations based on the SDE and ODE pharmacokinetic models are compared to real DCE-MRI data. They show that the use of SDE provides robustness in the estimation results. A simulation study confirms these results.     

Evolutionary optimization with data collocation for reverse engineering of biological networks

MOTIVATION: Modern experimental biology is moving away from analyses of single elements to whole-organism measurements. Such measured time-course data contain a wealth of information about the structure and dynamic of the pathway or network. The dynamic modeling of the whole systems is formulated as a reverse problem that requires a well-suited mathematical model and a very efficient computational method to identify the model structure and parameters. Numerical integration for differential equations and finding global parameter values are still two major challenges in this field of the parameter estimation of nonlinear dynamic biological systems. RESULTS: We compare three techniques of parameter estimation for nonlinear dynamic biological systems. In the proposed scheme, the modified collocation method is applied to convert the differential equations to the system of algebraic equations. The observed time-course data are then substituted into the algebraic system equations to decouple system interactions in order to obtain the approximate model profiles. Hybrid differential evolution (HDE) with population size of five is able to find a global solution. The method is not only suited for parameter estimation but also can be applied for structure identification. The solution obtained by HDE is then used as the starting point for a local search method to yield the refined estimates.     

Bayesian inference for dynamic transcriptional regulation; the Hes1 system as a case study

MOTIVATION: In this study, we address the problem of estimating the parameters of regulatory networks and provide the first application of Markov chain Monte Carlo (MCMC) methods to experimental data. As a case study, we consider a stochastic model of the Hes1 system expressed in terms of stochastic differential equations (SDEs) to which rigorous likelihood methods of inference can be applied. When fitting continuous-time stochastic models to discretely observed time series the lengths of the sampling intervals are important, and much of our study addresses the problem when the data are sparse. RESULTS: We estimate the parameters of an autoregulatory network providing results both for simulated and real experimental data from the Hes1 system. We develop an estimation algorithm using MCMC techniques which are flexible enough to allow for the imputation of latent data on a finer time scale and the presence of prior information about parameters which may be informed from other experiments as well as additional measurement error.     

Kinetic modeling of rhamnolipid production by Pseudomonas aeruginosa PAO1 including cell density-dependent regulation

The production of rhamnolipid biosurfactants by Pseudomonas aeruginosa is under complex control of a quorum sensing-dependent regulatory network. Due to a lack of understanding of the kinetics applicable to the process and relevant interrelations of variables, current processes for rhamnolipid production are based on heuristic approaches. To systematically establish a knowledge-based process for rhamnolipid production, a deeper understanding of the time-course and coupling of process variables is required. By combining reaction kinetics, stoichiometry, and experimental data, a process model for rhamnolipid production with P. aeruginosa PAO1 on sunflower oil was developed as a system of coupled ordinary differential equations (ODEs). In addition, cell density-based quorum sensing dynamics were included in the model. The model comprises a total of 36 parameters, 14 of which are yield coefficients and 7 of which are substrate affinity and inhibition constants. Of all 36 parameters, 30 were derived from dedicated experimental results, literature, and databases and 6 of them were used as fitting parameters. The model is able to describe data on biomass growth, substrates, and products obtained from a reference batch process and other validation scenarios. The model presented describes the time-course and interrelation of biomass, relevant substrates, and products on a process level while including a kinetic representation of cell density-dependent regulatory mechanisms.     

Parameter estimation of kinetic models from metabolic profiles: two-phase dynamic decoupling method

MOTIVATION: Time-series measurements of metabolite concentration have become increasingly more common, providing data for building kinetic models of metabolic networks using ordinary differential equations (ODEs). In practice, however, such time-course data are usually incomplete and noisy, and the estimation of kinetic parameters from these data is challenging. Practical limitations due to data and computational aspects, such as solving stiff ODEs and finding global optimal solution to the estimation problem, give motivations to develop a new estimation procedure that can circumvent some of these constraints. RESULTS: In this work, an incremental and iterative parameter estimation method is proposed that combines and iterates between two estimation phases. One phase involves a decoupling method, in which a subset of model parameters that are associated with measured metabolites, are estimated using the minimization of slope errors. Another phase follows, in which the ODE model is solved one equation at a time and the remaining model parameters are obtained by minimizing concentration errors. The performance of this two-phase method was tested on a generic branched metabolic pathway and the glycolytic pathway of Lactococcus lactis. The results showed that the method is efficient in getting accurate parameter estimates, even when some information is missing.     

Forcing Function Diagnostics for Nonlinear Dynamics

Summary .  This article investigates the problem of model diagnostics for systems described by nonlinear ordinary differential equations (ODEs). I propose modeling lack of fit as a time-varying correction to the right-hand side of a proposed differential equation. This correction can be described as being a set of additive forcing functions, estimated from data. Representing lack of fit in this manner allows us to graphically investigate model inadequacies and to suggest model improvements. I derive lack-of-fit tests based on estimated forcing functions. Model building in partially observed systems of ODEs is particularly difficult and I consider the problem of identification of forcing functions in these systems. The methods are illustrated with examples from computational neuroscience.     

Bifurcation analysis informs Bayesian inference in the Hes1 feedback loop

Background  

Ordinary differential equations (ODEs) are an important tool for describing the dynamics of biological systems. However, for ODE models to be useful, their parameters must first be calibrated. Parameter estimation, that is, finding parameter values given experimental data, is an inference problem that can be treated systematically through a Bayesian framework.     

Estimating parameters for generalized mass action models using constraint propagation

As modern molecular biology moves towards the analysis of biological systems as opposed to their individual components, the need for appropriate mathematical and computational techniques for understanding the dynamics and structure of such systems is becoming more pressing. For example, the modeling of biochemical systems using ordinary differential equations (ODEs) based on high-throughput, time-dense profiles is becoming more common-place, which is necessitating the development of improved techniques to estimate model parameters from such data. Due to the high dimensionality of this estimation problem, straight-forward optimization strategies rarely produce correct parameter values, and hence current methods tend to utilize genetic/evolutionary algorithms to perform non-linear parameter fitting. Here, we describe a completely deterministic approach, which is based on interval analysis. This allows us to examine entire sets of parameters, and thus to exhaust the global search within a finite number of steps. In particular, we show how our method may be applied to a generic class of ODEs used for modeling biochemical systems called Generalized Mass Action Models (GMAs). In addition, we show that for GMAs our method is amenable to the technique in interval arithmetic called constraint propagation, which allows great improvement of its efficiency. To illustrate the applicability of our method we apply it to some networks of biochemical reactions appearing in the literature, showing in particular that, in addition to estimating system parameters in the absence of noise, our method may also be used to recover the topology of these networks.     

A model of sequential branching in hierarchical cell fate determination

Multipotent stem or progenitor cells undergo a sequential series of binary fate decisions, which ultimately generate the diversity of differentiated cells. Efforts to understand cell fate control have focused on simple gene regulatory circuits that predict the presence of multiple stable states, bifurcations and switch-like transitions. However, existing gene network models do not explain more complex properties of cell fate dynamics such as the hierarchical branching of developmental paths. Here, we construct a generic minimal model of the genetic regulatory network controlling cell fate determination, which exhibits five elementary characteristics of cell differentiation: stability, directionality, branching, exclusivity, and promiscuous expression. We argue that a modular architecture comprising repeated network elements reproduces these features of differentiation by sequentially repressing selected modules and hence restricting the dynamics to lower dimensional subspaces of the high-dimensional state space. We implement our model both with ordinary differential equations (ODEs), to explore the role of bifurcations in producing the one-way character of differentiation, and with stochastic differential equations (SDEs), to demonstrate the effect of noise on the system. We further argue that binary cell fate decisions are prevalent in cell differentiation due to general features of the underlying dynamical system. This minimal model makes testable predictions about the structural basis for directional, discrete and diversifying cell phenotype development and thus can guide the evaluation of real gene regulatory networks that govern differentiation.     

Inference in dynamic systems using B‐splines and quasilinearized ODE penalties

Nonlinear (systems of) ordinary differential equations (ODEs) are common tools in the analysis of complex one‐dimensional dynamic systems. We propose a smoothing approach regularized by a quasilinearized ODE‐based penalty. Within the quasilinearized spline‐based framework, the estimation reduces to a conditionally linear problem for the optimization of the spline coefficients. Furthermore, standard ODE compliance parameter(s) selection criteria are applicable. We evaluate the performances of the proposed strategy through simulated and real data examples. Simulation studies suggest that the proposed procedure ensures more accurate estimates than standard nonlinear least squares approaches when the state (initial and/or boundary) conditions are not known.     

Numerical Simulation for the Unsteady MHD Flow and Heat Transfer of Couple Stress Fluid over a Rotating Disk

The present work is devoted to study the numerical simulation for unsteady MHD flow and heat transfer of a couple stress fluid over a rotating disk. A similarity transformation is employed to reduce the time dependent system of nonlinear partial differential equations (PDEs) to ordinary differential equations (ODEs). The Runge-Kutta method and shooting technique are employed for finding the numerical solution of the governing system. The influences of governing parameters viz. unsteadiness parameter, couple stress and various physical parameters on velocity, temperature and pressure profiles are analyzed graphically and discussed in detail.     

Efficient solution of ordinary differential equations modeling electrical activity in cardiac cells

The contraction of the heart is preceded and caused by a cellular electro-chemical reaction, causing an electrical field to be generated. Performing realistic computer simulations of this process involves solving a set of partial differential equations, as well as a large number of ordinary differential equations (ODEs) characterizing the reactive behavior of the cardiac tissue. Experiments have shown that the solution of the ODEs contribute significantly to the total work of a simulation, and there is thus a strong need to utilize efficient solution methods for this part of the problem. This paper presents how an efficient implicit Runge-Kutta method may be adapted to solve a complicated cardiac cell model consisting of 31 ODEs, and how this solver may be coupled to a set of PDE solvers to provide complete simulations of the electrical activity.     

Integrating genealogical and dynamical modelling to infer escape and reversion rates in HIV epitopes

The rates of escape and reversion in response to selection pressure arising from the host immune system, notably the cytotoxic T-lymphocyte (CTL) response, are key factors determining the evolution of HIV. Existing methods for estimating these parameters from cross-sectional population data using ordinary differential equations (ODEs) ignore information about the genealogy of sampled HIV sequences, which has the potential to cause systematic bias and overestimate certainty. Here, we describe an integrated approach, validated through extensive simulations, which combines genealogical inference and epidemiological modelling, to estimate rates of CTL escape and reversion in HIV epitopes. We show that there is substantial uncertainty about rates of viral escape and reversion from cross-sectional data, which arises from the inherent stochasticity in the evolutionary process. By application to empirical data, we find that point estimates of rates from a previously published ODE model and the integrated approach presented here are often similar, but can also differ several-fold depending on the structure of the genealogy. The model-based approach we apply provides a framework for the statistical analysis and hypothesis testing of escape and reversion in population data and highlights the need for longitudinal and denser cross-sectional sampling to enable accurate estimate of these key parameters.     

Analysis and Simulation of Division- and Label-Structured Population Models

In most biological studies and processes, cell proliferation and population dynamics play an essential role. Due to this ubiquity, a multitude of mathematical models has been developed to describe these processes. While the simplest models only consider the size of the overall populations, others take division numbers and labeling of the cells into account. In this work, we present a modeling and computational framework for proliferating cell populations undergoing symmetric cell division, which incorporates both the discrete division number and continuous label dynamics. Thus, it allows for the consideration of division number-dependent parameters as well as the direct comparison of the model prediction with labeling experiments, e.g., performed with Carboxyfluorescein succinimidyl ester (CFSE), and can be shown to be a generalization of most existing models used to describe these data. We prove that under mild assumptions the resulting system of coupled partial differential equations (PDEs) can be decomposed into a system of ordinary differential equations (ODEs) and a set of decoupled PDEs, which drastically reduces the computational effort for simulating the model. Furthermore, the PDEs are solved analytically and the ODE system is truncated, which allows for the prediction of the label distribution of complex systems using a low-dimensional system of ODEs. In addition to modeling the label dynamics, we link the label-induced fluorescence to the measure fluorescence which includes autofluorescence. Furthermore, we provide an analytical approximation for the resulting numerically challenging convolution integral. This is illustrated by modeling and simulating a proliferating population with division number-dependent proliferation rate.