Modeling and Estimation of Kinetic Parameters and Replicative Fitness of HIV-1 from Flow-Cytometry-Based Growth Competition Experiments

Growth competition assays have been developed to quantify the relative fitness of HIV-1 mutants. In this article, we develop mathematical models to describe viral/cellular dynamic interactions in the assay system from which the competitive fitness indices or parameters are defined. In our previous HIV-viral fitness experiments, the concentration of uninfected target cells was assumed to be constant (Wu et al. 2006). But this may not be true in some experiments. In addition, dual infection may frequently occur in viral fitness experiments and may not be ignorable. Here, we relax these two assumptions and extend our earlier viral fitness model (Wu et al. 2006). The resulting models then become nonlinear ODE systems for which closed-form solutions are not achievable. In the new model, the viral relative fitness is a function of time since it depends on the target cell concentration. First, we studied the structure identifiability of the nonlinear ODE models. The identifiability analysis showed that all parameters in the proposed models are identifiable from the flow-cytometry-based experimental data that we collected. We then employed a global optimization approach (the differential evolution algorithm) to directly estimate the kinetic parameters as well as the relative fitness index in the nonlinear ODE models using nonlinear least square regression based on the experimental data. Practical identifiability was investigated via Monte Carlo simulations.     

Hierarchical Bayesian inference for HIV dynamic differential equation models incorporating multiple treatment factors

Studies on HIV dynamics in AIDS research are very important in understanding the pathogenesis of HIV‐1 infection and also in assessing the effectiveness of antiretroviral (ARV) treatment. Viral dynamic models can be formulated through a system of nonlinear ordinary differential equations (ODE), but there has been only limited development of statistical methodologies for inference. This article, motivated by an AIDS clinical study, discusses a hierarchical Bayesian nonlinear mixed‐effects modeling approach to dynamic ODE models without a closed‐form solution. In this model, we fully integrate viral load, medication adherence, drug resistance, pharmacokinetics, baseline covariates and time‐dependent drug efficacy into the data analysis for characterizing long‐term virologic responses. Our method is implemented by a data set from an AIDS clinical study. The results suggest that modeling HIV dynamics and virologic responses with consideration of time‐varying clinical factors as well as baseline characteristics may be important for HIV/AIDS studies in providing quantitative guidance to better understand the virologic responses to ARV treatment and to help the evaluation of clinical trial design in response to existing therapies.     

A nonlinear,second-order population model

Delay differential, difference, and partial differential equation models are being used more extensively to explain single-species population oscillations and limit cycle behavior. Ordinary differential equation (ODE) models have been largely ignored. This is because first-order ODE models are inherently monotonic. Certainly this is not usual population behavior in the real world. If it is assumed that the per capita growth rate of a population changes over time as a result of regulating factors impinging on it, then a more realistic biological model results. The model translates into a second-order nonlinear ODE. Such a model can exhibit oscillatory and limit cycle as well as monotonic solutions, i.e., behavior for which non-ODE models have been used to explain. Although first-order ODE models are gross simplifications of real phenomena, ODE models in general should not be disregarded as important analytical tools.     

Maximum likelihood estimation in dynamical models of HIV

The study of dynamical models of HIV infection, based on a system of nonlinear ordinary differential equations (ODE), has considerably improved the knowledge of its pathogenesis. While the first models used simplified ODE systems and analyzed each patient separately, recent works dealt with inference in non-simplified models borrowing strength from the whole sample. The complexity of these models leads to great difficulties for inference and only the Bayesian approach has been attempted by now. We propose a full likelihood inference, adapting a Newton-like algorithm for these particular models. We consider a relatively complex ODE model for HIV infection and a model for the observations including the issue of detection limits. We apply this approach to the analysis of a clinical trial of antiretroviral therapy (ALBI ANRS 070) and we show that the whole algorithm works well in a simulation study.     

Quantifying Immune Response to Influenza Virus Infection via Multivariate Nonlinear ODE Models with Partially Observed State Variables and Time-Varying Parameters

Influenza A virus (IAV) infection continues to be a global health threat, as evidenced by the outbreak of the novel A/California/7/2009 IAV strain. Previous flu vaccines have proven less effective than hoped for emerging IAV strains, indicating a more thorough understanding of immune responses to primary infection is needed. One issue is the difficulty in directly measuring many key parameters and variables of the immune response. To address these issues, we considered a comprehensive workflow for statistical inference for ordinary differential question (ODE) models with partially observed variables and time-varying parameters, including identifiability analysis, two-stage and NLS estimation, model selection, etc. In particular, we proposed a novel one-step method to verify parameter identifiability and formulate estimating equations simultaneously. Thus, the pseudo-LS method can now deal with general ODE models with partially observed state variables for the first time. Using this workflow, we verified the relative significance of various immune factors to virus control, including target epithelial cells, cytotoxic T-lymphocyte (CD8+) cells, and IAV specific antibodies (IgG and IgM). Factors other than cytotoxic T-lymphocyte (CTL) killing contributed the most to the loss of infected epithelial cells, though the effects of CTL are still significant. IgM antibody was found to be the major contributor to neutralization of free infectious viral particles. Also, the maximum viral load, which correlates well with mortality, was found to depend more on viral replication rates than infectivity. In contrast to current hypotheses, the results obtained via our methods suggest that IgM antibody and viral replication rates may be worth of further explorations in vaccine development.     

A comparison of statistical methods for detecting differentially expressed genes from RNA-seq data

RNA-Seq technologies are quickly revolutionizing genomic studies, and statistical methods for RNA-seq data are under continuous development. Timely review and comparison of the most recently proposed statistical methods will provide a useful guide for choosing among them for data analysis. Particular interest surrounds the ability to detect differential expression (DE) in genes. Here we compare four recently proposed statistical methods, edgeR, DESeq, baySeq, and a method with a two-stage Poisson model (TSPM), through a variety of simulations that were based on different distribution models or real data. We compared the ability of these methods to detect DE genes in terms of the significance ranking of genes and false discovery rate control. All methods compared are implemented in freely available software. We also discuss the availability and functions of the currently available versions of these software.     

Backward bifurcation and oscillations in a nested immuno-eco-epidemiological model

This paper introduces a novel partial differential equation immuno-eco-epidemiological model of competition in which one species is affected by a disease while another can compete with it directly and by lowering the first species' immune response to the infection, a mode of competition termed stress-induced competition. When the disease is chronic, and the within-host dynamics are rapid, we reduce the partial differential equation model (PDE) to a three-dimensional ordinary differential equation (ODE) model. The ODE model exhibits backward bifurcation and sustained oscillations caused by the stress-induced competition. Furthermore, the ODE model, although not a special case of the PDE model, is useful for detecting backward bifurcation and oscillations in the PDE model. Backward bifurcation related to stress-induced competition allows the second species to persist for values of its invasion number below one. Furthermore, stress-induced competition leads to destabilization of the coexistence equilibrium and sustained oscillations in the PDE model. We suggest that complex systems such as this one may be studied by appropriately designed simple ODE models.     

A unique transformation from ordinary differential equations to reaction networks

Many models in Systems Biology are described as a system of Ordinary Differential Equations, which allows for transient, steady-state or bifurcation analysis when kinetic information is available. Complementary structure-related qualitative analysis techniques have become increasingly popular in recent years, like qualitative model checking or pathway analysis (elementary modes, invariants, flux balance analysis, graph-based analyses, chemical organization theory, etc.). They do not rely on kinetic information but require a well-defined structure as stochastic analysis techniques equally do. In this article, we look into the structure inference problem for a model described by a system of Ordinary Differential Equations and provide conditions for the uniqueness of its solution. We describe a method to extract a structured reaction network model, represented as a bipartite multigraph, for example, a continuous Petri net (CPN), from a system of Ordinary Differential Equations (ODEs). A CPN uniquely defines an ODE, and each ODE can be transformed into a CPN. However, it is not obvious under which conditions the transformation of an ODE into a CPN is unique, that is, when a given ODE defines exactly one CPN. We provide biochemically relevant sufficient conditions under which the derived structure is unique and counterexamples showing the necessity of each condition. Our method is implemented and available; we illustrate it on some signal transduction models from the BioModels database. A prototype implementation of the method is made available to modellers at http://contraintes.inria.fr/~soliman/ode2pn.html, and the data mentioned in the "Results" section at http://contraintes.inria.fr/~soliman/ode2pn_data/. Our results yield a new recommendation for the import/export feature of tools supporting the SBML exchange format.     

Numerical discretization-based estimation methods for ordinary differential equation models via penalized spline smoothing with applications in biomedical research

Differential equations are extensively used for modeling dynamics of physical processes in many scientific fields such as engineering, physics, and biomedical sciences. Parameter estimation of differential equation models is a challenging problem because of high computational cost and high-dimensional parameter space. In this article, we propose a novel class of methods for estimating parameters in ordinary differential equation (ODE) models, which is motivated by HIV dynamics modeling. The new methods exploit the form of numerical discretization algorithms for an ODE solver to formulate estimating equations. First, a penalized-spline approach is employed to estimate the state variables and the estimated state variables are then plugged in a discretization formula of an ODE solver to obtain the ODE parameter estimates via a regression approach. We consider three different order of discretization methods, Euler's method, trapezoidal rule, and Runge-Kutta method. A higher-order numerical algorithm reduces numerical error in the approximation of the derivative, which produces a more accurate estimate, but its computational cost is higher. To balance the computational cost and estimation accuracy, we demonstrate, via simulation studies, that the trapezoidal discretization-based estimate is the best and is recommended for practical use. The asymptotic properties for the proposed numerical discretization-based estimators are established. Comparisons between the proposed methods and existing methods show a clear benefit of the proposed methods in regards to the trade-off between computational cost and estimation accuracy. We apply the proposed methods t an HIV study to further illustrate the usefulness of the proposed approaches.     

Simple testing procedures for the Holling type II model

In this paper, we consider predator–prey data that can be viewed as solutions to a planar system of ordinary differential equations (ODE) observed with random error. The ODE system admits a limit cycle, while the random error is supposed to act additively in the log-scale. One of the oldest such systems is Holling’s type II model. In spite of its simplicity, it is still very popular in data analyses, although more sophisticated models have been introduced in the literature. We propose a simple way of deciding whether a set of predator–prey pairs is indicative or not of a departure from this basic model by exploiting the geometric properties of the solution in the phase plane. To illustrate our method, we use simulated and real data.     

Stochastic simulations on a model of circadian rhythm generation

Biological phenomena are often modeled by differential equations, where states of a model system are described by continuous real values. When we consider concentrations of molecules as dynamical variables for a set of biochemical reactions, we implicitly assume that numbers of the molecules are large enough so that their changes can be regarded as continuous and they are described deterministically. However, for a system with small numbers of molecules, changes in their numbers are apparently discrete and molecular noises become significant. In such cases, models with deterministic differential equations may be inappropriate, and the reactions must be described by stochastic equations. In this study, we focus a clock gene expression for a circadian rhythm generation, which is known as a system involving small numbers of molecules. Thus it is appropriate for the system to be modeled by stochastic equations and analyzed by methodologies of stochastic simulations. The interlocked feedback model proposed by Ueda et al. as a set of deterministic ordinary differential equations provides a basis of our analyses. We apply two stochastic simulation methods, namely Gillespie's direct method and the stochastic differential equation method also by Gillespie, to the interlocked feedback model. To this end, we first reformulated the original differential equations back to elementary chemical reactions. With those reactions, we simulate and analyze the dynamics of the model using two methods in order to compare them with the dynamics obtained from the original deterministic model and to characterize dynamics how they depend on the simulation methodologies.     

PLMaddon: a power-law module for the Matlab SBToolbox

PLMaddon is a General Public License (GPL) software module designed to expand the current version of the SBToolbox (a Matlab toolbox for systems biology; www.sbtoolbox.org) with a set of functions for the analysis of power-law models, a specific class of kinetic models, set in ordinary differential equations (ODE) and in which the kinetic orders can have positive/negative non-integer values. The module includes functions to generate power-law Taylor expansions of other ODE models (e.g. Michaelis-Menten type models), as well as algorithms to estimate steady-states. The robustness and sensitivity of the models can also be analysed and visualized by computing the power-law's logarithmic gains and sensitivities.     

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.     

Establishing Informative Prior for Gene Expression Variance from Public Databases

Identifying differential expressed genes across various conditions or genotypes is the most typical approach to studying the regulation of gene expression. An estimate of gene-specific variance is often needed for the assessment of statistical significance in most differential expression (DE) detection methods, including linear models (e.g., for transformed and normalized microarray data) and generalized linear models (e.g., for count data in RNAseq). Due to a common limit in sample size, the variance estimate is often unstable in small experiments. Shrinkage estimates using empirical Bayes methods have proven useful in improving the variance estimate, hence improving the detection of DE. The most widely used empirical Bayes methods borrow information across genes within the same experiments. In these methods, genes are considered exchangeable or exchangeable conditioning on expression level. We propose, with the increasing accumulation of expression data, borrowing information from historical data on the same gene can provide better estimate of gene-specific variance, thus further improve DE detection. Specifically, we show that the variation of gene expression is truly gene-specific and reproducible between different experiments. We present a new method to establish informative gene-specific prior on the variance of expression using existing public data, and illustrate how to shrink the variance estimate and detect DE. We demonstrate improvement in DE detection under our strategy compared to leading DE detection methods.     

A dynamic model for functional mapping of biological rhythms

Functional mapping is a statistical method for mapping quantitative trait loci (QTLs) that regulate the dynamic pattern of a biological trait. This method integrates mathematical aspects of biological complexity into a mixture model for genetic mapping and tests the genetic effects of QTLs by comparing genotype-specific curve parameters. As a way of quantitatively specifying the dynamic behaviour of a system, differential equations have proved to be powerful for modelling and unravelling the biochemical, molecular, and cellular mechanisms of a biological process, such as biological rhythms. The equipment of functional mapping with biologically meaningful differential equations provides new insights into the genetic control of any dynamic processes. We formulate a new functional mapping framework for a dynamic biological rhythm by incorporating a group of ordinary differential equations (ODE). The Runge–Kutta fourth-order algorithm was implemented to estimate the parameters that define the system of ODE. The new model will find its implications for understanding the interplay between gene interactions and developmental pathways in complex biological rhythms.     

A dynamic model for functional mapping of biological rhythms

Functional mapping is a statistical method for mapping quantitative trait loci (QTLs) that regulate the dynamic pattern of a biological trait. This method integrates mathematical aspects of biological complexity into a mixture model for genetic mapping and tests the genetic effects of QTLs by comparing genotype-specific curve parameters. As a way of quantitatively specifying the dynamic behavior of a system, differential equations have proven to be powerful for modeling and unraveling the biochemical, molecular, and cellular mechanisms of a biological process, such as biological rhythms. The equipment of functional mapping with biologically meaningful differential equations provides new insights into the genetic control of any dynamic processes. We formulate a new functional mapping framework for a dynamic biological rhythm by incorporating a group of ordinary differential equations (ODE). The Runge-Kutta fourth order algorithm was implemented to estimate the parameters that define the system of ODE. The new model will find its implications for understanding the interplay between gene interactions and developmental pathways in complex biological rhythms.     

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.