Parameter identifiability and estimation of HIV/AIDS dynamic models

A novel within-host model of HIV and nutrition

In this paper we develop a four compartment within-host model of nutrition and HIV. We show that the model has two equilibria: an infection-free equilibrium and infection equilibrium. The infection free equilibrium is locally asymptotically stable when the basic reproduction number $ \mathcal{R}_0 < 1 $, and unstable when $ \mathcal{R}_0 > 1 $. The infection equilibrium is locally asymptotically stable if $ \mathcal{R}_0 > 1 $ and an additional condition holds. We show that the within-host model of HIV and nutrition is structured to reveal its parameters from the observations of viral load, CD4 cell count and total protein data. We then estimate the model parameters for these 3 data sets. We have also studied the practical identifiability of the model parameters by performing Monte Carlo simulations, and found that the rate of clearance of the virus by immunoglobulins is practically unidentifiable, and that the rest of the model parameters are only weakly identifiable given the experimental data. Furthermore, we have studied how the data frequency impacts the practical identifiability of model parameters.

How robust are estimates of key parameters in standard viral dynamic models?

Mathematical models of viral infection have been developed, fitted to data, and provide insight into disease pathogenesis for multiple agents that cause chronic infection, including HIV, hepatitis C, and B virus. However, for agents that cause acute infections or during the acute stage of agents that cause chronic infections, viral load data are often collected after symptoms develop, usually around or after the peak viral load. Consequently, we frequently lack data in the initial phase of viral growth, i.e., when pre-symptomatic transmission events occur. Missing data may make estimating the time of infection, the infectious period, and parameters in viral dynamic models, such as the cell infection rate, difficult. However, having extra information, such as the average time to peak viral load, may improve the robustness of the estimation. Here, we evaluated the robustness of estimates of key model parameters when viral load data prior to the viral load peak is missing, when we know the values of some parameters and/or the time from infection to peak viral load. Although estimates of the time of infection are sensitive to the quality and amount of available data, particularly pre-peak, other parameters important in understanding disease pathogenesis, such as the loss rate of infected cells, are less sensitive. Viral infectivity and the viral production rate are key parameters affecting the robustness of data fits. Fixing their values to literature values can help estimate the remaining model parameters when pre-peak data is missing or limited. We find a lack of data in the pre-peak growth phase underestimates the time to peak viral load by several days, leading to a shorter predicted growth phase. On the other hand, knowing the time of infection (e.g., from epidemiological data) and fixing it results in good estimates of dynamical parameters even in the absence of early data. While we provide ways to approximate model parameters in the absence of early viral load data, our results also suggest that these data, when available, are needed to estimate model parameters more precisely.

Minimally sufficient experimental design using identifiability analysis

Mathematical models are increasingly being developed and calibrated in tandem with data collection, empowering scientists to intervene in real time based on quantitative model predictions. Well-designed experiments can help augment the predictive power of a mathematical model but the question of when to collect data to maximize its utility for a model is non-trivial. Here we define data as model-informative if it results in a unique parametrization, assessed through the lens of practical identifiability. The framework we propose identifies an optimal experimental design (how much data to collect and when to collect it) that ensures parameter identifiability (permitting confidence in model predictions), while minimizing experimental time and costs. We demonstrate the power of the method by applying it to a modified version of a classic site-of-action pharmacokinetic/pharmacodynamic model that describes distribution of a drug into the tumor microenvironment (TME), where its efficacy is dependent on the level of target occupancy in the TME. In this context, we identify a minimal set of time points when data needs to be collected that robustly ensures practical identifiability of model parameters. The proposed methodology can be applied broadly to any mathematical model, allowing for the identification of a minimally sufficient experimental design that collects the most informative data.

A computational approach to identifiability analysis for a model of the propagation and control of COVID-19 in Chile

A computational approach is adapted to analyze the parameter identifiability of a compartmental model. The model is intended to describe the progression of the COVID-19 pandemic in Chile during the initial phase in early 2020 when government declared quarantine measures. The computational approach to analyze the structural and practical identifiability is applied in two parts, one for synthetic data and another for some Chilean regional data. The first part defines the identifiable parameter sets when these recover the true parameters used to create the synthetic data. The second part compares the results derived from synthetic data, estimating the identifiable parameter sets from regional Chilean epidemic data. Experiments provide evidence of the loss of identifiability if some initial conditions are estimated, the period of time used to fit is before the peak, and if a significant proportion of the population is involved in quarantine periods.

A Continuation Technique for Maximum Likelihood Estimators in Biological Models

Estimating model parameters is a crucial step in mathematical modelling and typically involves minimizing the disagreement between model predictions and experimental data. This calibration data can change throughout a study, particularly if modelling is performed simultaneously with the calibration experiments, or during an on-going public health crisis as in the case of the COVID-19 pandemic. Consequently, the optimal parameter set, or maximal likelihood estimator (MLE), is a function of the experimental data set. Here, we develop a numerical technique to predict the evolution of the MLE as a function of the experimental data. We show that, when considering perturbations from an initial data set, our approach is significantly more computationally efficient that re-fitting model parameters while producing acceptable model fits to the updated data. We use the continuation technique to develop an explicit functional relationship between fit model parameters and experimental data that can be used to measure the sensitivity of the MLE to experimental data. We then leverage this technique to select between model fits with similar information criteria, a priori determine the experimental measurements to which the MLE is most sensitive, and suggest additional experiment measurements that can resolve parameter uncertainty.

HIV infection dynamics and viral rebound: Modeling results from humanized mice

Despite years of combined antiretroviral therapy (cART), HIV persists in infected individuals. The virus also rebounds after the cessation of cART. The sources contributing to viral persistence and rebound are not fully understood. When viral rebound occurs, what affects the time to rebound and how to delay the rebound remain unclear. In this paper, we started with the data fitting of an HIV infection model to the viral load data in treated and untreated humanized myeloid-only mice (MoM) in which macrophages serve as the target of HIV infection. By fixing the parameter values for macrophages from the MoM fitting, we fit a mathematical model including the infection of two target cell populations to the viral load data from humanized bone marrow/liver/thymus (BLT) mice, in which both CD4+ T cells and macrophages are the target of HIV infection. Data fitting suggests that the viral load decay in BLT mice under treatment has three phases. The loss of infected CD4+ T cells and macrophages is a major contributor to the first two phases of viral decay, and the last phase may be due to the latent infection of CD4+ T cells. Numerical simulations using parameter estimates from the data fitting show that the pre-ART viral load and the latent reservoir size at treatment cessation can affect viral growth rate and predict the time to viral rebound. Model simulations further reveal that early and prolonged cART can delay the viral rebound after cessation of treatment, which may have implications in the search for functional control of HIV infection.

From data to interpretable models: machine learning for soil moisture forecasting

Soil moisture is critical to agricultural business, ecosystem health, and certain hydrologically driven natural disasters. Monitoring data, though, is prone to instrumental noise, wide ranging extrema, and nonstationary response to rainfall where ground conditions change. Furthermore, existing soil moisture models generally forecast poorly for time periods greater than a few hours. To improve such forecasts, we introduce two data-driven models, the Naive Accumulative Representation (NAR) and the Additive Exponential Accumulative Representation (AEAR). Both of these models are rooted in deterministic, physically based hydrology, and we study their capabilities in forecasting soil moisture over time periods longer than a few hours. Learned model parameters represent the physically based unsaturated hydrological redistribution processes of gravity and suction. We validate our models using soil moisture and rainfall time series data collected from a steep gradient, post-wildfire site in southern California. Data analysis is complicated by rapid landscape change observed in steep, burned hillslopes in response to even small to moderate rain events. The proposed NAR and AEAR models are, in forecasting experiments, shown to be competitive with several established and state-of-the-art baselines. The AEAR model fits the data well for three distinct soil textures at variable depths below the ground surface (5, 15, and 30 cm). Similar robust results are demonstrated in controlled, laboratory-based experiments. Our AEAR model includes readily interpretable hydrologic parameters and provides more accurate forecasts than existing models for time horizons of 10-24 h. Such extended periods of warning for natural disasters, such as floods and landslides, provide actionable knowledge to reduce loss of life and property.

Predicting and forecasting the impact of local outbreaks of COVID-19: use of SEIR-D quantitative epidemiological modelling for healthcare demand and capacity

BACKGROUND

The world is experiencing local/regional hotspots and spikes in the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which causes COVID-19 disease. We aimed to formulate an applicable epidemiological model to accurately predict and forecast the impact of local outbreaks of COVID-19 to guide the local healthcare demand and capacity, policy-making and public health decisions.

METHODS

The model utilized the aggregated daily COVID-19 situation reports (including counts of daily admissions, discharges and bed occupancy) from the local National Health Service (NHS) hospitals and COVID-19-related weekly deaths in hospitals and other settings in Sussex (population 1.7 million), Southeast England. These data sets corresponded to the first wave of COVID-19 infections from 24 March to 15 June 2020. A novel epidemiological predictive and forecasting model was then derived based on the local/regional surveillance data. Through a rigorous inverse parameter inference approach, the model parameters were estimated by fitting the model to the data in an optimal sense and then subsequent validation.

RESULTS

The inferred parameters were physically reasonable and matched up to the widely used parameter values derived from the national data sets by Biggerstaff M, Cowling BJ, Cucunubá ZM et al. (Early insights from statistical and mathematical modeling of key epidemiologic parameters of COVID-19, Emerging infectious diseases. 2020;26(11)). We validate the predictive power of our model by using a subset of the available data and comparing the model predictions for the next 10, 20 and 30 days. The model exhibits a high accuracy in the prediction, even when using only as few as 20 data points for the fitting.

CONCLUSIONS

We have demonstrated that by using local/regional data, our predictive and forecasting model can be utilized to guide the local healthcare demand and capacity, policy-making and public health decisions to mitigate the impact of COVID-19 on the local population. Understanding how future COVID-19 spikes/waves could possibly affect the regional populations empowers us to ensure the timely commissioning and organization of services. The flexibility of timings in the model, in combination with other early-warning systems, produces a time frame for these services to prepare and isolate capacity for likely and potential demand within regional hospitals. The model also allows local authorities to plan potential mortuary capacity and understand the burden on crematoria and burial services. The model algorithms have been integrated into a web-based multi-institutional toolkit, which can be used by NHS hospitals, local authorities and public health departments in other regions of the UK and elsewhere. The parameters, which are locally informed, form the basis of predicting and forecasting exercises accounting for different scenarios and impacts of COVID-19 transmission.

Testing structural identifiability by a simple scaling method

Successful mathematical modeling of biological processes relies on the expertise of the modeler to capture the essential mechanisms in the process at hand and on the ability to extract useful information from empirical data. A model is said to be structurally unidentifiable, if different quantitative sets of parameters provide the same observable outcome. This is typical (but not exclusive) of partially observed problems in which only a few variables can be experimentally measured. Most of the available methods to test the structural identifiability of a model are either too complex mathematically for the general practitioner to be applied, or require involved calculations or numerical computation for complex non-linear models. In this work, we present a new analytical method to test structural identifiability of models based on ordinary differential equations, based on the invariance of the equations under the scaling transformation of its parameters. The method is based on rigorous mathematical results but it is easy and quick to apply, even to test the identifiability of sophisticated highly non-linear models. We illustrate our method by example and compare its performance with other existing methods in the literature.

Parameter and State Estimation in a Cholera Model with Threshold Immunology: A Case Study of Senegal

It is often impossible to measure all states affecting spread of a disease. In cholera, asymptomatic and cholera pathogen densities are not practically measurable despite playing a big role in its transmission. They are referred to as inaccessible states of the model and can only be manipulated using the measurable states of the given model. Our interest lies in estimating such states and the parameters catalyzing the spread. A mathematical model for cholera dynamics consisting of five compartments (susceptible, symptomatic, asymptomatic, recovered and bacteria population) with a minimum infection dose (MID) is considered. A method based on observer (from modern control theory) is proposed to estimate the state variables not accessible to measurement and the time-dependent parameters from real data of Senegal. We suppose that the total population of Senegal, the monthly reported cholera-induced deaths and the monthly recovered individuals are known inputs obtainable from real data, and the monthly new cholera cases the system output. An auxiliary system is used, an observer whose solutions converge exponentially to those of an original system and solely utilize known inputs and output of the model. Thus, the estimation of the unmeasured state variables like the pathogen density and the yearly asymptomatic population within the human community playing an important role in the spread of cholera is possible. We derive the expressions for time-dependent infection rate, induced cholera death rate and symptomatic recovery rate and their estimations done using real data. Numerical simulations are then performed for the validation of estimation results. The system together with the observer designed is detectable but is not observable. The observer delivered estimates reflect a close trend already ascertained by other researchers. They indicate the existence of bacteria and asymptomatic individuals in the population of Senegal at any single time for the duration of collection of the data. The ever existence of cholera pathogen explains the endemicity of Cholera in Senegal and other sub-Saharan-African countries owing to role played by the asymptomatic individual in the bacteria density. As such, the heath authorities in Senegal need to educate the general public on hygiene irrespective of observable symptoms to lower the possible number of new infections. We have analytically showed and numerically confirmed the exponential convergence to zero of the estimation errors resulting from the observer model hence the high quality of the estimates.

Modelling Degradation and Replication Kinetics of the Zika Virus In Vitro Infection

Mathematical models of in vitro viral kinetics help us understand and quantify the main determinants underlying the virus–host cell interactions. We aimed to provide a numerical characterization of the Zika virus (ZIKV) in vitro infection kinetics, an arthropod-borne emerging virus that has gained public recognition due to its association with microcephaly in newborns. The mathematical model of in vitro viral infection typically assumes that degradation of extracellular infectious virus proceeds in an exponential manner, that is, each viral particle has the same probability of losing infectivity at any given time. We incubated ZIKV stock in the cell culture media and sampled with high frequency for quantification over the course of 96 h. The data showed a delay in the virus degradation in the first 24 h followed by a decline, which could not be captured by the model with exponentially distributed decay time of infectious virus. Thus, we proposed a model, in which inactivation of infectious ZIKV is gamma distributed and fit the model to the temporal measurements of infectious virus remaining in the media. The model was able to reproduce the data well and yielded the decay time of infectious ZIKV to be 40 h. We studied the in vitro ZIKV infection kinetics by conducting cell infection at two distinct multiplicity of infection and measuring viral loads over time. We fit the mathematical model of in vitro viral infection with gamma distributed degradation time of infectious virus to the viral growth data and identified the timespans and rates involved within the ZIKV-host cell interplay. Our mathematical analysis combined with the data provides a well-described example of non-exponential viral decay dynamics and presents numerical characterization of in vitro infection with ZIKV.

Model Averaging in Viral Dynamic Models

The paucity of experimental data makes both inference and prediction particularly challenging in viral dynamic models. In the presence of several candidate models, a common strategy is model selection (MS), in which models are fitted to the data but only results obtained with the "best model" are presented. However, this approach ignores model uncertainty, which may lead to inaccurate predictions. When several models provide a good fit to the data, another approach is model averaging (MA) that weights the predictions of each model according to its consistency to the data. Here, we evaluated by simulations in a nonlinear mixed-effect model framework the performances of MS and MA in two realistic cases of acute viral infection, i.e., (1) inference in the presence of poorly identifiable parameters, namely, initial viral inoculum and eclipse phase duration, (2) uncertainty on the mechanisms of action of the immune response. MS was associated in some scenarios with a large rate of false selection. This led to a coverage rate lower than the nominal coverage rate of 0.95 in the majority of cases and below 0.50 in some scenarios. In contrast, MA provided better estimation of parameter uncertainty, with coverage rates ranging from 0.72 to 0.98 and mostly comprised within the nominal coverage rate. Finally, MA provided similar predictions than those obtained with MS. In conclusion, parameter estimates obtained with MS should be taken with caution, especially when several models well describe the data. In this situation, MA has better performances and could be performed to account for model uncertainty.

Herpes simplex virus-2 dynamics as a probe to measure the extremely rapid and spatially localized tissue-resident T-cell response

Herpes simplex virus-2 infection is characterized by frequent episodic shedding in the genital tract. Expansion in HSV-2 viral load early during episodes is extremely rapid. However, the virus invariably peaks within 18 hours and is eliminated nearly as quickly. A critical feature of HSV-2 shedding episodes is their heterogeneity. Some episodes peak at 10⁸ HSV DNA copies, last for weeks due to frequent viral re-expansion, and lead to painful ulcers, while others only reach 10³ HSV DNA copies and are eliminated within hours and without symptoms. Within single micro-environments of infection, tissue-resident CD8+ T cells (T_RM ) appear to contain infection within a few days. Here, we review components of T_RM biology relevant to immune surveillance between HSV-2 shedding episodes and containment of infection upon detection of HSV-2 cognate antigen. We then describe the use of mathematical models to correlate large spatial gradients in T_RM density with the heterogeneity of observed shedding within a single person. We describe how models have been leveraged for clinical trial simulation, as well as future plans to model the interactions of multiple cellular subtypes within mucosa, predict the mechanism of action of therapeutic vaccines, and describe the dynamics of 3-dimensional infection environment during the natural evolution of an HSV-2 lesion.

Parameter estimation and mathematical modeling for the quantitative description of therapy failure due to drug resistance in gastrointestinal stromal tumor metastasis to the liver

In this work we develop a general mathematical model and devise a practical identifiability approach for gastrointestinal stromal tumor (GIST) metastasis to the liver, with the aim of quantitatively describing therapy failure due to drug resistance. To this end, we have modeled metastatic growth and therapy failure produced by resistance to two standard treatments based on tyrosine kinase inhibitors (Imatinib and Sunitinib) that have been observed clinically in patients with GIST metastasis to the liver. The parameter identification problem is difficult to solve, since there are no general results on this issue for models based on ordinary differential equations (ODE) like the ones studied here. We propose a general modeling framework based on ODE for GIST metastatic growth and therapy failure due to drug resistance and analyzed five different model variants, using medical image observations (CT scans) from patients that exhibit drug resistance. The associated parameter estimation problem was solved using the Nelder-Mead simplex algorithm, by adding a regularization term to the objective function to address model instability, and assessing the agreement of either an absolute or proportional error in the objective function. We compared the goodness of fit to data for the proposed model variants, as well as evaluated both error forms in order to improve parameter estimation results. From the model variants analyzed, we identified the one that provides the best fit to all the available patient data sets, as well as the best assumption in computing the objective function (absolute or proportional error). This is the first work that reports mathematical models capable of capturing and quantitatively describing therapy failure due to drug resistance based on clinical images in a patient-specific manner.

Correlation Between Anti-gp41 Antibodies and Virus Infectivity Decay During Primary HIV-1 Infection

Recent experiments have suggested that the infectivity of simian immunodeficiency virus (SIV) and human immunodeficiency virus type-1 (HIV-1) in plasma decreases over time during primary infection. Because anti-gp41 antibodies are produced early during HIV-1 infection and form antibody-virion complexes, we studied if such early HIV-1 specific antibodies are correlated with the decay in HIV-1 infectivity. Using a viral dynamic model that allows viral infectivity to decay and frequent early viral load data obtained from 6 plasma donors we estimate that HIV-1 infectivity begins to decay after about 2 weeks of infection. The length of this delay is consistent with the time before antibody-virion complexes were detected in the plasma of these donors and is correlated (p = 0.023, r = 0.87) with the time for antibodies to be first detected in plasma. Importantly, we identify that the rate of infectivity decay is significantly correlated with the rate of increase in plasma anti-gp41 IgG concentration (p = 0.046, r = 0.82) and the increase in IgM+IgG anti-gp41 concentration (p = 8.37 × 10⁻⁴, r = 0.98). Furthermore, we found that the viral load decay after the peak did not have any significant correlation with the rate of anti-gp41 IgM or IgG increase. These results indicate that early anti-gp41 antibodies may cause viral infectivity decay, but may not contribute significantly to controlling post-peak viral load, likely due to insufficient quantity or affinity. Our findings may be helpful to devise strategies, including antibody-based vaccines, to control acute HIV-1 infection.

Dynamical compensation and structural identifiability of biological models: Analysis, implications, and reconciliation

The concept of dynamical compensation has been recently introduced to describe the ability of a biological system to keep its output dynamics unchanged in the face of varying parameters. However, the original definition of dynamical compensation amounts to lack of structural identifiability. This is relevant if model parameters need to be estimated, as is often the case in biological modelling. Care should we taken when using an unidentifiable model to extract biological insight: the estimated values of structurally unidentifiable parameters are meaningless, and model predictions about unmeasured state variables can be wrong. Taking this into account, we explore alternative definitions of dynamical compensation that do not necessarily imply structural unidentifiability. Accordingly, we show different ways in which a model can be made identifiable while exhibiting dynamical compensation. Our analyses enable the use of the new concept of dynamical compensation in the context of parameter identification, and reconcile it with the desirable property of structural identifiability.

A robust behaviour is a desirable feature in many biological systems. The study of mechanisms capable of maintaining the transient response unchanged despite environmental disturbances has recently motivated the introduction of a new concept: Dynamical Compensation (DC). However, the original definition of DC with respect to a parameter amounts to structural unidentifiability of that parameter, which means that it cannot be estimated by measuring the model output. Since most biological models have unknown parameters that need to be estimated, DC can be considered a negative property for the purpose of model identification. In this paper we reconcile these two conflicting views by proposing a new definition of DC that captures its intended biological meaning (i.e. robustness, which should be a systemic property, intrinsic to the dynamics) while making it distinct from structural unidentifiability (which is a modelling property that depends on decisions made by the modeller, such as the choice of model outputs or unknown parameters, and on experimental constraints). Our definition enables a model to have DC with respect to a structurally identifiable parameter, thus increasing the applicability of the concept.

Quantifying rotavirus kinetics in the REH tumor cell line using in vitro data

Globally, rotavirus is the most common cause of diarrhea in children younger than 5 years of age, however, a quantitative understanding of the infection dynamics is still lacking. In this paper, we present the first study to extract viral kinetic parameters for in vitro rotavirus infections in the REH cell tumor line. We use a mathematical model of viral kinetics to extract parameter values by fitting the model to data from rotavirus infection of REH cells. While accurate results for some of the parameters of the mathematical model were not achievable due to its global non-identifiability, we are able to quantify approximately the time course of the infection for the first time. We also find that the basic reproductive number of rotavirus, which gives the number of secondary infections from a single infected cell, is much greater than one. Quantifying the kinetics of rotavirus leads not only to a better understanding of the infection process, but also provides a method for quantitative comparison of kinetics of different strains or for quantifying the effectiveness of antiviral treatment.

Dynamics of lentiviral infection in vivo in the absence of adaptive immune responses

Understanding the dynamics of acute viral infection is crucial for developing strategies to prevent and control infection. In this study, lentiviral dynamics in a host without adaptive immunity were examined in order to determine kinetic parameters of infection and quantify the effect of neutralizing antibodies in preventing infection, using mathematical modeling of data from equine infectious anemia virus (EIAV) infection of horses with severe combined immunodeficiency (SCID). Estimated parameters were used to calculate the basic reproductive number and virus doubling time and found that the rate that antibodies neutralized virus was ~18 times greater than the virus clearance rate. These results establish EIAV replication kinetics in SCID horses and the minimal efficacy of antibodies that blocked infection. Furthermore, they indicate that EIAV is at most mildly cytopathic. This study advances our understanding of EIAV infection and may have important implications for the control of other viral infections, including HIV.

Study design and parameter estimability for spatial and temporal ecological models

The statistical tools available to ecologists are becoming increasingly sophisticated, allowing more complex, mechanistic models to be fit to ecological data. Such models have the potential to provide new insights into the processes underlying ecological patterns, but the inferences made are limited by the information in the data. Statistical nonestimability of model parameters due to insufficient information in the data is a problem too‐often ignored by ecologists employing complex models. Here, we show how a new statistical computing method called data cloning can be used to inform study design by assessing the estimability of parameters under different spatial and temporal scales of sampling. A case study of parasite transmission from farmed to wild salmon highlights that assessing the estimability of ecologically relevant parameters should be a key step when designing studies in which fitting complex mechanistic models is the end goal.

Analysis of Practical Identifiability of a Viral Infection Model

Mathematical modelling approaches have granted a significant contribution to life sciences and beyond to understand experimental results. However, incomplete and inadequate assessments in parameter estimation practices hamper the parameter reliability, and consequently the insights that ultimately could arise from a mathematical model. To keep the diligent works in modelling biological systems from being mistrusted, potential sources of error must be acknowledged. Employing a popular mathematical model in viral infection research, existing means and practices in parameter estimation are exemplified. Numerical results show that poor experimental data is a main source that can lead to erroneous parameter estimates despite the use of innovative parameter estimation algorithms. Arbitrary choices of initial conditions as well as data asynchrony distort the parameter estimates but are often overlooked in modelling studies. This work stresses the existence of several sources of error buried in reports of modelling biological systems, voicing the need for assessing the sources of error, consolidating efforts in solving the immediate difficulties, and possibly reconsidering the use of mathematical modelling to quantify experimental data.

Parameter and state estimation in a Neisseria meningitidis model: A study case of Niger

Neisseria meningitidis (Nm) is a major cause of bacterial meningitidis outbreaks in Africa and the Middle East. The availability of yearly reported meningitis cases in the African meningitis belt offers the opportunity to analyze the transmission dynamics and the impact of control strategies. In this paper, we propose a method for the estimation of state variables that are not accessible to measurements and an unknown parameter in a Nm model. We suppose that the yearly number of Nm induced mortality and the total population are known inputs, which can be obtained from data, and the yearly number of new Nm cases is the model output. We also suppose that the Nm transmission rate is an unknown parameter. We first show how the recruitment rate into the population can be estimated using real data of the total population and Nm induced mortality. Then, we use an auxiliary system called observer whose solutions converge exponentially to those of the original model. This observer does not use the unknown infection transmission rate but only uses the known inputs and the model output. This allows us to estimate unmeasured state variables such as the number of carriers that play an important role in the transmission of the infection and the total number of infected individuals within a human community. Finally, we also provide a simple method to estimate the unknown Nm transmission rate. In order to validate the estimation results, numerical simulations are conducted using real data of Niger.

Hierarchical effects of pro-inflammatory cytokines on the post-influenza susceptibility to pneumococcal coinfection

In the course of influenza A virus (IAV) infections, a secondary bacterial infection frequently leads to serious respiratory conditions provoking high hospitalization and death tolls. Although abundant pro-inflammatory responses have been reported as key contributing factors for these severe dual infections, the relative contributions of cytokines remain largely unclear. In the current study, mathematical modelling based on murine experimental data dissects IFN-γ as a cytokine candidate responsible for impaired bacterial clearance, thereby promoting bacterial growth and systemic dissemination during acute IAV infection. We also found a time-dependent detrimental role of IL-6 in curtailing bacterial outgrowth which was not as distinct as for IFN-γ. Our numerical simulations suggested a detrimental effect of IFN-γ alone and in synergism with IL-6 but no conclusive pathogenic effect of IL-6 and TNF-α alone. This work provides a rationale to understand the potential impact of how to manipulate temporal immune components, facilitating the formulation of hypotheses about potential therapeutic strategies to treat coinfections.

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.

A model-based initial guess for estimating parameters in systems of ordinary differential equations

The inverse problem of parameter estimation from noisy observations is a major challenge in statistical inference for dynamical systems. Parameter estimation is usually carried out by optimizing some criterion function over the parameter space. Unless the optimization process starts with a good initial guess, the estimation may take an unreasonable amount of time, and may converge to local solutions, if at all. In this article, we introduce a novel technique for generating good initial guesses that can be used by any estimation method. We focus on the fairly general and often applied class of systems linear in the parameters. The new methodology bypasses numerical integration and can handle partially observed systems. We illustrate the performance of the method using simulations and apply it to real data.

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, and 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.

Ebola virus infection modeling and identifiability problems

The recent outbreaks of Ebola virus (EBOV) infections have underlined the impact of the virus as a major threat for human health. Due to the high biosafety classification of EBOV (level 4), basic research is very limited. Therefore, the development of new avenues of thinking to advance quantitative comprehension of the virus and its interaction with the host cells is urgently needed to tackle this lethal disease. Mathematical modeling of the EBOV dynamics can be instrumental to interpret Ebola infection kinetics on quantitative grounds. To the best of our knowledge, a mathematical modeling approach to unravel the interaction between EBOV and the host cells is still missing. In this paper, a mathematical model based on differential equations is used to represent the basic interactions between EBOV and wild-type Vero cells in vitro. Parameter sets that represent infectivity of pathogens are estimated for EBOV infection and compared with influenza virus infection kinetics. The average infecting time of wild-type Vero cells by EBOV is slower than in influenza infection. Simulation results suggest that the slow infecting time of EBOV could be compensated by its efficient replication. This study reveals several identifiability problems and what kind of experiments are necessary to advance the quantification of EBOV infection. A first mathematical approach of EBOV dynamics and the estimation of standard parameters in viral infections kinetics is the key contribution of this work, paving the way for future modeling works on EBOV infection.

Generalized Ordinary Differential Equation Models

Existing estimation methods for ordinary differential equation (ODE) models are not applicable to discrete data. The generalized ODE (GODE) model is therefore proposed and investigated for the first time. We develop the likelihood-based parameter estimation and inference methods for GODE models. We propose robust computing algorithms and rigorously investigate the asymptotic properties of the proposed estimator by considering both measurement errors and numerical errors in solving ODEs. The simulation study and application of our methods to an influenza viral dynamics study suggest that the proposed methods have a superior performance in terms of accuracy over the existing ODE model estimation approach and the extended smoothing-based (ESB) method.

Estimation of nonlinear differential equation model for glucose-insulin dynamics in type I diabetic patients using generalized smoothing

In this work, we develop an ordinary differential equations (ODE) model of physiological regulation of glycemia in type 1 diabetes mellitus (T1DM) patients in response to meals and intravenous insulin infusion. Unlike for majority of existing mathematical models of glucose-insulin dynamics, parameters in our model are estimable from a relatively small number of noisy observations of plasma glucose and insulin concentrations. For estimation, we adopt the generalized smoothing estimation of nonlinear dynamic systems of Ramsay et al. (2007). In this framework, the ODE solution is approximated with a penalized spline, where the ODE model is incorporated in the penalty. We propose to optimize the generalized smoothing by using penalty weights that minimize the covariance penalties criterion (Efron, 2004). The covariance penalties criterion provides an estimate of the prediction error for nonlinear estimation rules resulting from nonlinear and/or non-homogeneous ODE models, such as our model of glucose-insulin dynamics. We also propose to select the optimal number and location of knots for B-spline bases used to represent the ODE solution. The results of the small simulation study demonstrate advantages of optimized generalized smoothing in terms of smaller estimation errors for ODE parameters and smaller prediction errors for solutions of differential equations. Using the proposed approach to analyze the glucose and insulin concentration data in T1DM patients we obtained good approximation of global glucose-insulin dynamics and physiologically meaningful parameter estimates.

Identification and experimental validation of an HIV model for HAART treated patients

The objective of this paper is to identify the parameters of a human immunodeficiency virus (HIV) evolution model from a clinical data set of patients treated with two different highly active antiretroviral therapy (HAART) protocols. After introducing a model with six state variables, a preliminary step considers the reduction of the number of parameters to be identified by means of sensitivity analysis, and then identifiability items are discussed. A nonlinear optimization-based procedure for identification is developed, which divides the unknown parameters into two families: the group dependent and the patient dependent parameters. Numerical results show that the identified model can be individually adapted to each patient and this result is promising for predicting the effects (e.g., failures or successes) of therapeutic actions.

Influence of a priori Information, Designs, and Undetectable Data on Individual Parameters Estimation and Prediction of Hepatitis C Treatment Outcome

Hepatitis C viral kinetic analysis based on nonlinear mixed effect models can be used to individualize treatment. For that purpose, it is necessary to obtain precise estimation of individual parameters. Here, we evaluated by simulation the influence on Bayesian individual parameter estimation and outcome prediction of a priori information on population parameters, viral load sampling designs, and methods for handling data below detection limit (BDL). We found that a precise estimation of both individual parameters and treatment outcome could be obtained using as few as six measurements in the first month of therapy. This result remained valid even when incorrect a priori information on population parameters was set as long as the parameters were identifiable and BDL data were properly handled. However, setting wrong values for a priori population parameters could lead to severe estimation/prediction errors if BDL data were ignored and not properly accounted in the likelihood function.

A practical approach to parameter estimation applied to model predicting heart rate regulation

Mathematical models have long been used for prediction of dynamics in biological systems. Recently, several efforts have been made to render these models patient specific. One way to do so is to employ techniques to estimate parameters that enable model based prediction of observed quantities. Knowledge of variation in parameters within and between groups of subjects have potential to provide insight into biological function. Often it is not possible to estimate all parameters in a given model, in particular if the model is complex and the data is sparse. However, it may be possible to estimate a subset of model parameters reducing the complexity of the problem. In this study, we compare three methods that allow identification of parameter subsets that can be estimated given a model and a set of data. These methods will be used to estimate patient specific parameters in a model predicting baroreceptor feedback regulation of heart rate during head-up tilt. The three methods include: structured analysis of the correlation matrix, analysis via singular value decomposition followed by QR factorization, and identification of the subspace closest to the one spanned by eigenvectors of the model Hessian. Results showed that all three methods facilitate identification of a parameter subset. The "best" subset was obtained using the structured correlation method, though this method was also the most computationally intensive. Subsets obtained using the other two methods were easier to compute, but analysis revealed that the final subsets contained correlated parameters. In conclusion, to avoid lengthy computations, these three methods may be combined for efficient identification of parameter subsets.

Modelling HIV-1 2-LTR dynamics following raltegravir intensification

A model of reservoir activation and viral replication is introduced accounting for the production of 2-LTR HIV-1 DNA circles following antiviral intensification with the HIV integrase inhibitor raltegravir, considering contributions of de novo infection events and exogenous sources of infected cells, including quiescent infected cell activation. The model shows that a monotonic increase in measured 2-LTR concentration post intensification is consistent with limited de novo infection primarily maintained by sources of infected cells unaffected by raltegravir, such as quiescent cell activation, while a transient increase in measured 2-LTR concentration is consistent with significant levels of efficient (R0 > 1) de novo infection. The model is validated against patient data from the INTEGRAL study and is shown to have a statistically significant fit relative to the null hypothesis of random measurement variation about a mean. We obtain estimates and confidence intervals for the model parameters, including 2-LTR half-life. Seven of the 13 patients with detectable 2-LTR concentrations from the INTEGRAL study have measured 2-LTR dynamics consistent with significant levels of efficient replication of the virus prior to treatment intensification.

Modelling the course of an HIV infection: insights from ecology and evolution

The Human Immunodeficiency Virus (HIV) is one of the most threatening viral agents. This virus infects approximately 33 million people, many of whom are unaware of their status because, except for flu-like symptoms right at the beginning of the infection during the acute phase, the disease progresses more or less symptom-free for 5 to 10 years. During this asymptomatic phase, the virus slowly destroys the immune system until the onset of AIDS when opportunistic infections like pneumonia or Kaposi’s sarcoma can overcome immune defenses. Mathematical models have played a decisive role in estimating important parameters (e.g., virion clearance rate or life-span of infected cells). However, most models only account for the acute and asymptomatic latency phase and cannot explain the progression to AIDS. Models that account for the whole course of the infection rely on different hypotheses to explain the progression to AIDS. The aim of this study is to review these models, present their technical approaches and discuss the robustness of their biological hypotheses. Among the few models capturing all three phases of an HIV infection, we can distinguish between those that mainly rely on population dynamics and those that involve virus evolution. Overall, the modeling quest to capture the dynamics of an HIV infection has improved our understanding of the progression to AIDS but, more generally, it has also led to the insight that population dynamics and evolutionary processes can be necessary to explain the course of an infection.

Which of our modeling predictions are robust?

In theoretical ecology it is well known that the steady state expressions of the variables in a food chain crucially depend on the parity of the length of the chain. This poses a major problem for modeling real food webs because it is difficult to establish their true number of trophic levels, with sometimes rare predators and often rampant pathogens. Similar problems arise in the modeling of chronic viral infections. We review examples where seemingly general interpretations strongly depend on the number of levels in a model, and on its specific equations. This Perspective aims to open the discussion on this problem.

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.

Parameter estimation and uncertainty quantification for an epidemic model

We examine estimation of the parameters of Susceptible-Infective-Recovered (SIR) models in the context of least squares. We review the use of asymptotic statistical theory and sensitivity analysis to obtain measures of uncertainty for estimates of the model parameters and the basic reproductive number (R0)---an epidemiologically significant parameter grouping. We find that estimates of different parameters, such as the transmission parameter and recovery rate, are correlated, with the magnitude and sign of this correlation depending on the value of R0. Situations are highlighted in which this correlation allows R0 to be estimated with greater ease than its constituent parameters. Implications of correlation for parameter identifiability are discussed. Uncertainty estimates and sensitivity analysis are used to investigate how the frequency at which data is sampled affects the estimation process and how the accuracy and uncertainty of estimates improves as data is collected over the course of an outbreak. We assess the informativeness of individual data points in a given time series to determine when more frequent sampling (if possible) would prove to be most beneficial to the estimation process. This technique can be used to design data sampling schemes in more general contexts.

Quantitative analysis of viral persistence and transient viral load rebound from HIV clinical data

Highly active antiretroviral therapy (HAART) suppresses HIV RNA viral load below the limit of detection for many patients. However, clinical data demonstrates that the HIV virus is not eradicated by HAART, even in patients who have had no detectable virus for 7 years [1]. One possible reason is that a stable resting latent reservoir with a long half-life exists in resting memory CD4(+)T cells [2]. In this paper, we propose a mathematical model with a constant contribution of a stable latent reservoir and identified this constant by using one patient's data from AutoVac HAART interruption study [3]. Many patients also have transient rebounds of plasma viral RNA (viral blips) under otherwise successful control of the virus by HAART. Activation of latently infected cells can explain these transient rebounds of viral load. Little quantitative analysis about the activation of reservoir has been done based on any clinical experiment data. Here, we model the activation dynamics of the reservoir by a time-independent activation rate and estimate this rate by using the clinical data from the AutoVac HAART interruption study [3].

Parameter Estimation for Differential Equation Models Using a Framework of Measurement Error in Regression Models

Differential equation (DE) models are widely used in many scientific fields that include engineering, physics and biomedical sciences. The so-called "forward problem", the problem of simulations and predictions of state variables for given parameter values in the DE models, has been extensively studied by mathematicians, physicists, engineers and other scientists. However, the "inverse problem", the problem of parameter estimation based on the measurements of output variables, has not been well explored using modern statistical methods, although some least squares-based approaches have been proposed and studied. In this paper, we propose parameter estimation methods for ordinary differential equation models (ODE) based on the local smoothing approach and a pseudo-least squares (PsLS) principle under a framework of measurement error in regression models. The asymptotic properties of the proposed PsLS estimator are established. We also compare the PsLS method to the corresponding SIMEX method and evaluate their finite sample performances via simulation studies. We illustrate the proposed approach using an application example from an HIV dynamic study.

Structural identifiability of systems biology models: a critical comparison of methods

Analysing the properties of a biological system through in silico experimentation requires a satisfactory mathematical representation of the system including accurate values of the model parameters. Fortunately, modern experimental techniques allow obtaining time-series data of appropriate quality which may then be used to estimate unknown parameters. However, in many cases, a subset of those parameters may not be uniquely estimated, independently of the experimental data available or the numerical techniques used for estimation. This lack of identifiability is related to the structure of the model, i.e. the system dynamics plus the observation function. Despite the interest in knowing a priori whether there is any chance of uniquely estimating all model unknown parameters, the structural identifiability analysis for general non-linear dynamic models is still an open question. There is no method amenable to every model, thus at some point we have to face the selection of one of the possibilities. This work presents a critical comparison of the currently available techniques. To this end, we perform the structural identifiability analysis of a collection of biological models. The results reveal that the generating series approach, in combination with identifiability tableaus, offers the most advantageous compromise among range of applicability, computational complexity and information provided.

Joint modeling of the clinical progression and of the biomarkers' dynamics using a mechanistic model

Joint models are used to rigorously explore the relationship between the dynamics of biomarkers and clinical events. In the context of HIV infection, where the multivariate dynamics of HIV-RNA and CD4 are complex, a mechanistic approach based on a system of nonlinear differential equations naturally takes into account the correlation between the biomarkers. Using data from a randomized clinical trial comparing dual antiretroviral therapy to a single drug regimen, a full maximum likelihood approach is proposed to explore the relationship between the evolution of the biomarkers and the time to a clinical event. The role of each marker as an independent predictor of disease progression is assessed. We show that the joint dynamics of HIV-RNA and CD4 captures the effect of antiretroviral treatment; the CD4 dynamics alone is found to capture most but not all of the treatment effect.

Maximum likelihood estimation of long-term HIV dynamic models and antiviral response

HIV dynamics studies, based on differential equations, have significantly improved the knowledge on HIV infection. While first studies used simplified short-term dynamic models, recent works considered more complex long-term models combined with a global analysis of whole patient data based on nonlinear mixed models, increasing the accuracy of the HIV dynamic analysis. However statistical issues remain, given the complexity of the problem. We proposed to use the SAEM (stochastic approximation expectation-maximization) algorithm, a powerful maximum likelihood estimation algorithm, to analyze simultaneously the HIV viral load decrease and the CD4 increase in patients using a long-term HIV dynamic system. We applied the proposed methodology to the prospective COPHAR2-ANRS 111 trial. Very satisfactory results were obtained with a model with latent CD4 cells defined with five differential equations. One parameter was fixed, the 10 remaining parameters (eight with between-patient variability) of this model were well estimated. We showed that the efficacy of nelfinavir was reduced compared to indinavir and lopinavir.

Assessing the in vitro fitness of an oseltamivir-resistant seasonal A/H1N1 influenza strain using a mathematical model

In 2007, the A/Brisbane/59/2007 (H1N1) seasonal influenza virus strain acquired the oseltamivir-resistance mutation H275Y in its neuraminidase (NA) gene. Although previous studies had demonstrated that this mutation impaired the replication capacity of the influenza virus in vitro and in vivo, the A/Brisbane/59/2007 H275Y oseltamivir-resistant mutant completely out-competed the wild-type (WT) strain and was, in the 2008–2009 influenza season, the primary A/H1N1 circulating strain. Using a combination of plaque and viral yield assays, and a simple mathematical model, approximate values were extracted for two basic viral kinetics parameters of the in vitro infection. In the ST6GalI-MDCK cell line, the latent infection period (i.e., the time for a newly infected cell to start releasing virions) was found to be 1–3 h for the WT strain and more than 7 h for the H275Y mutant. The infecting time (i.e., the time for a single infectious cell to cause the infection of another one) was between 30 and 80 min for the WT, and less than 5 min for the H275Y mutant. Single-cycle viral yield experiments have provided qualitative confirmation of these findings. These results, though preliminary, suggest that the increased fitness success of the A/Brisbane/59/2007 H275Y mutant may be due to increased infectivity compensating for an impaired or delayed viral release, and are consistent with recent evidence for the mechanistic origins of fitness reduction and recovery in NA expression. The method applied here can reconcile seemingly contradictory results from the plaque and yield assays as two complementary views of replication kinetics, with both required to fully capture a strain's fitness.