Differential equation modeling of HIV viral fitness experiments: model identification, model selection, and multimodel inference

Summary .  Many biological processes and systems can be described by a set of differential equation (DE) models. However, literature in statistical inference for DE models is very sparse. We propose statistical estimation, model selection, and multimodel averaging methods for HIV viral fitness experiments in vitro that can be described by a set of nonlinear ordinary differential equations (ODE). The parameter identifiability of the ODE models is also addressed. We apply the proposed methods and techniques to experimental data of viral fitness for HIV-1 mutant 103N. We expect that the proposed modeling and inference approaches for the DE models can be widely used for a variety of biomedical studies.     

Estimating a predator-prey dynamical model with the parameter cascades method

Summary .   Ordinary differential equations (ODEs) are widely used in ecology to describe the dynamical behavior of systems of interacting populations. However, systems of ODEs rarely provide quantitative solutions that are close to real field observations or experimental data because natural systems are subject to environmental and demographic noise and ecologists are often uncertain about the correct parameterization. In this article we introduce "parameter cascades" as an improved method to estimate ODE parameters such that the corresponding ODE solutions fit the real data well. This method is based on the modified penalized smoothing with the penalty defined by ODEs and a generalization of profiled estimation, which leads to fast estimation and good precision for ODE parameters from noisy data. This method is applied to a set of ODEs originally developed to describe an experimental predator–prey system that undergoes oscillatory dynamics. The new parameterization considerably improves the fit of the ODE model to the experimental data sets. At the same time, our method reveals that important structural assumptions that underlie the original ODE model are essentially correct. The mathematical formulations of the two nonlinear interaction terms (functional responses) that link the ODEs in the predator–prey model are validated by estimating the functional responses nonparametrically from the real data. We suggest two major applications of "parameter cascades" to ecological modeling: It can be used to estimate parameters when original data are noisy, missing, or when no reliable priori estimates are available; it can help to validate the structural soundness of the mathematical modeling approach.     

Parameter estimation for stiff equations of biosystems using radial basis function networks

Background  

The modeling of dynamic systems requires estimating kinetic parameters from experimentally measured time-courses. Conventional global optimization methods used for parameter estimation, e.g. genetic algorithms (GA), consume enormous computational time because they require iterative numerical integrations for differential equations. When the target model is stiff, the computational time for reaching a solution increases further.     

HIV model parameter estimates from interruption trial data including drug efficacy and reservoir dynamics

Mathematical models based on ordinary differential equations (ODE) have had significant impact on understanding HIV disease dynamics and optimizing patient treatment. A model that characterizes the essential disease dynamics can be used for prediction only if the model parameters are identifiable from clinical data. Most previous parameter identification studies for HIV have used sparsely sampled data from the decay phase following the introduction of therapy. In this paper, model parameters are identified from frequently sampled viral-load data taken from ten patients enrolled in the previously published AutoVac HAART interruption study, providing between 69 and 114 viral load measurements from 3-5 phases of viral decay and rebound for each patient. This dataset is considerably larger than those used in previously published parameter estimation studies. Furthermore, the measurements come from two separate experimental conditions, which allows for the direct estimation of drug efficacy and reservoir contribution rates, two parameters that cannot be identified from decay-phase data alone. A Markov-Chain Monte-Carlo method is used to estimate the model parameter values, with initial estimates obtained using nonlinear least-squares methods. The posterior distributions of the parameter estimates are reported and compared for all patients.     

Targeted maximum likelihood estimation of the parameter of a marginal structural model

Targeted maximum likelihood estimation is a versatile tool for estimating parameters in semiparametric and nonparametric models. We work through an example applying targeted maximum likelihood methodology to estimate the parameter of a marginal structural model. In the case we consider, we show how this can be easily done by clever use of standard statistical software. We point out differences between targeted maximum likelihood estimation and other approaches (including estimating function based methods). The application we consider is to estimate the effect of adherence to antiretroviral medications on virologic failure in HIV positive individuals.     

Estimating growth and mortality in stage-structured populations

This paper presents a practical numerical method for separatingand estimating growth and mortality coefficients in stage- orsize-structured populations using only observations of the relativeor absolute abundance of each stage. The method involves writinga system of linear ordinary differential equations (ODEs) modellingthe rate of change of abundance. The solution of the differentialsystem can be numerically approximated using standard (e.g.sixth-order Runge-Kutta-Felhberg) methods. An optimization problemwhose solutions yield ‘optimal’ coefficients fora given model is formulated. The ODE numerical integration techniquecan then be employed to furnish required function and gradientinformation to the optimization algorithm. The data-fittingsoftware package ODRPACK is then successfully employed to estimateoptimal coefficients for the ODE population model. Simulationexperiments with four- and eight-stage model populations illustratethat the method results in the successful estimation of coefficientsof mortality and growth from abundance data.     

Numerical methods and parameter estimation of a structured population model with discrete events in the life history

We consider two numerical methods for the solution of a physiologically structured population (PSP) model with multiple life stages and discrete event reproduction. The model describes the dynamic behaviour of a predator-prey system consisting of rotifers predating on algae. The nitrate limited algal prey population is modelled unstructured and described by an ordinary differential equation (ODE). The formulation of the rotifer dynamics is based on a simple physiological model for their two life stages, the egg and the adult stage. An egg is produced when an energy buffer reaches a threshold value. The governing equations are coupled partial differential equations (PDE) with initial and boundary conditions. The population models together with the equation for the dynamics of the nutrient result in a chemostat model. Experimental data are used to estimate the model parameters. The results obtained with the explicit finite difference (FD) technique compare well with those of the Escalator Boxcar Train (EBT) method. This justifies the use of the fast FD method for the parameter estimation, a procedure which involves repeated solution of the model equations.     

Novel recurrent neural network for modelling biological networks: Oscillatory p53 interaction dynamics

Understanding the control of cellular networks consisting of gene and protein interactions and their emergent properties is a central activity of Systems Biology research. For this, continuous, discrete, hybrid, and stochastic methods have been proposed. Currently, the most common approach to modelling accurate temporal dynamics of networks is ordinary differential equations (ODE). However, critical limitations of ODE models are difficulty in kinetic parameter estimation and numerical solution of a large number of equations, making them more suited to smaller systems. In this article, we introduce a novel recurrent artificial neural network (RNN) that addresses above limitations and produces a continuous model that easily estimates parameters from data, can handle a large number of molecular interactions and quantifies temporal dynamics and emergent systems properties. This RNN is based on a system of ODEs representing molecular interactions in a signalling network. Each neuron represents concentration change of one molecule represented by an ODE. Weights of the RNN correspond to kinetic parameters in the system and can be adjusted incrementally during network training. The method is applied to the p53-Mdm2 oscillation system – a crucial component of the DNA damage response pathways activated by a damage signal. Simulation results indicate that the proposed RNN can successfully represent the behaviour of the p53-Mdm2 oscillation system and solve the parameter estimation problem with high accuracy. Furthermore, we presented a modified form of the RNN that estimates parameters and captures systems dynamics from sparse data collected over relatively large time steps. We also investigate the robustness of the p53-Mdm2 system using the trained RNN under various levels of parameter perturbation to gain a greater understanding of the control of the p53-Mdm2 system. Its outcomes on robustness are consistent with the current biological knowledge of this system. As more quantitative data become available on individual proteins, the RNN would be able to refine parameter estimation and mapping of temporal dynamics of individual signalling molecules as well as signalling networks as a system. Moreover, RNN can be used to modularise large signalling networks.     

Transport Effects on Multiple-Component Reactions in Optical Biosensors

Optical biosensors are often used to measure kinetic rate constants associated with chemical reactions. Such instruments operate in the surface–volume configuration, in which ligand molecules are convected through a fluid-filled volume over a surface to which receptors are confined. Currently, scientists are using optical biosensors to measure the kinetic rate constants associated with DNA translesion synthesis—a process critical to DNA damage repair. Biosensor experiments to study this process involve multiple interacting components on the sensor surface. This multiple-component biosensor experiment is modeled with a set of nonlinear integrodifferential equations (IDEs). It is shown that in physically relevant asymptotic limits these equations reduce to a much simpler set of ordinary differential equations (ODEs). To verify the validity of our ODE approximation, a numerical method for the IDE system is developed and studied. Results from the ODE model agree with simulations of the IDE model, rendering our ODE model useful for parameter estimation.     

A novel cost function to estimate parameters of oscillatory biochemical systems

Oscillatory pathways are among the most important classes of biochemical systems with examples ranging from circadian rhythms and cell cycle maintenance. Mathematical modeling of these highly interconnected biochemical networks is needed to meet numerous objectives such as investigating, predicting and controlling the dynamics of these systems. Identifying the kinetic rate parameters is essential for fully modeling these and other biological processes. These kinetic parameters, however, are not usually available from measurements and most of them have to be estimated by parameter fitting techniques. One of the issues with estimating kinetic parameters in oscillatory systems is the irregularities in the least square (LS) cost function surface used to estimate these parameters, which is caused by the periodicity of the measurements. These irregularities result in numerous local minima, which limit the performance of even some of the most robust global optimization algorithms. We proposed a parameter estimation framework to address these issues that integrates temporal information with periodic information embedded in the measurements used to estimate these parameters. This periodic information is used to build a proposed cost function with better surface properties leading to fewer local minima and better performance of global optimization algorithms. We verified for three oscillatory biochemical systems that our proposed cost function results in an increased ability to estimate accurate kinetic parameters as compared to the traditional LS cost function. We combine this cost function with an improved noise removal approach that leverages periodic characteristics embedded in the measurements to effectively reduce noise. The results provide strong evidence on the efficacy of this noise removal approach over the previous commonly used wavelet hard-thresholding noise removal methods. This proposed optimization framework results in more accurate kinetic parameters that will eventually lead to biochemical models that are more precise, predictable, and controllable.     

Censored linear regression in the presence or absence of auxiliary survival information

There has been a rising interest in better exploiting auxiliary summary information from large databases in the analysis of smaller-scale studies that collect more comprehensive patient-level information. The purpose of this paper is twofold: first, we propose a novel approach to synthesize information from both the aggregate summary statistics and the individual-level data in censored linear regression. We show that the auxiliary information amounts to a system of nonsmooth estimating equations and thus can be combined with the conventional weighted log-rank estimating equations by using the generalized method of moments (GMM) approach. The proposed methodology can be further extended to account for the potential inconsistency in information from different sources. Second, in the absence of auxiliary information, we propose to improve estimation efficiency by combining the overidentified weighted log-rank estimating equations with different weight functions via the GMM framework. To deal with the nonsmooth GMM-type objective functions, we develop an asymptotics-guided algorithm for parameter and variance estimation. We establish the asymptotic normality of the proposed GMM-type estimators. Simulation studies show that the proposed estimators can yield substantial efficiency gain over the conventional weighted log-rank estimators. The proposed methods are applied to a pancreatic cancer study for illustration.     

Metal reduction kinetics in Shewanella

MOTIVATION: Metal reduction kinetics have been studied in cultures of dissimilatory metal reducing bacteria which include the Shewanella oneidensis strain MR-1. Estimation of system parameters from time-series data faces obstructions in the implementation depending on the choice of the mathematical model that captures the observed dynamics. The modeling of metal reduction is often based on Michaelis-Menten equations. These models are often developed using initial in vitro reaction rates and seldom match with in vivo reduction profiles. RESULTS: For metal reduction studies, we propose a model that is based on the power law representation that is effectively applied to the kinetics of metal reduction. The method yields reasonable parameter estimates and is illustrated with the analysis of time-series data that describes the dynamics of metal reduction in S.oneidensis strain MR-1. In addition, mixed metal studies involving the reduction of Uranyl (U(VI)) to the relatively insoluble tetravalent form (U(IV)) by S. alga strain (BR-Y) were studied in the presence of environmentally relevant iron hydrous oxides. For mixed metals, parameter estimation and curve fitting are accomplished with a generalized least squares formulation that handles systems of ordinary differential equations and is implemented in Matlab. It consists of an optimization algorithm (Levenberg-Marquardt, LSQCURVEFIT) and a numerical ODE solver. Simulation with the estimated parameters indicates that the model captures the experimental data quite well. The model uses the estimated parameters to predict the reduction rates of metals and mixed metals at varying concentrations. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.     

Median regression models for longitudinal data with dropouts

Summary .  Recently, median regression models have received increasing attention. When continuous responses follow a distribution that is quite different from a normal distribution, usual mean regression models may fail to produce efficient estimators whereas median regression models may perform satisfactorily. In this article, we discuss using median regression models to deal with longitudinal data with dropouts. Weighted estimating equations are proposed to estimate the median regression parameters for incomplete longitudinal data, where the weights are determined by modeling the dropout process. Consistency and the asymptotic distribution of the resultant estimators are established. The proposed method is used to analyze a longitudinal data set arising from a controlled trial of HIV disease ( Volberding et al., 1990 , The New England Journal of Medicine 322, 941–949). Simulation studies are conducted to assess the performance of the proposed method under various situations. An extension to estimation of the association parameters is outlined.     

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.     

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.     

A simple method for estimating the parameter of substitution rate variation among sites

When the rate variation among sites is described by a gamma distribution, an important problem is how to estimate the shape parameter alpha, which is an index of the degree of among-site rate variation. The parsimony-based methods for estimating alpha are simple but biased, i.e., alpha tends to be overestimated. On the other hand, the likelihood-based methods are asymptotically unbiased but take a huge amount of computational time. In this paper, we have developed a new method to solve this problem: we first estimate the expected number of substitutions at each site, which is corrected for multiple hits, and then estimate the parameter alpha. Our method is computationally as fast as the parsimony method, and the estimation accuracy is much higher than that of parsimony and similar to that of the likelihood method.      

Symbolic-numeric estimation of parameters in biochemical models by quantifier elimination

The sequencing of complete genomes allows analyses of the interactions between various biological molecules on a genomic scale, which prompted us to simulate the global behaviors of biological phenomena on the molecular level. One of the basic mathematical problems in the simulation is the parameter optimization in the kinetic model for complex dynamics, and many estimation methods have been designed. We introduce a new approach to estimate the parameters in biological kinetic models by quantifier elimination (QE), in combination with numerical simulation methods. The estimation method was applied to a model for the inhibition kinetics of HIV proteinase with ten parameters and nine variables, and attained the goodness of fit to 300 points of observed data with the same magnitude as that obtained by the previous estimation methods, remarkably by using only one or two points of data. Furthermore, the utilization of QE demonstrated the feasibility of the present method for elucidating the behavior of the parameters and the variables in the analyzed model. Therefore, the present symbolic-numeric method is a powerful approach to reveal the fundamental mechanisms of kinetic models, in addition to being a computational engine.