Robust Parameter Estimation of Rational Ordinary Differential Equations for NLME and QSP
Chris Rackauckas (1,2,3), Paul Lang (1), Elisabeth Roesch (1), Sebastian Micluta-Câmpeanu (1), Yosef Berman (4), Soo Go (4), Hoon Hong (4), Ilia Ilmer (4), Alexey Ovchinnikov (4), Pedro Soto (4), and Chee Yap (4)
JuliaHub (1), Pumas-AI (2), Computer Science and Artificial Intelligence Laboratory (CSAIL), Massachusetts Institute of Technology (3), New York University Courant Institute (4)
Introduction/Objectives:
Parameter estimation of ordinary differential equation (ODE) models is the essential step of many inference methods in nonlinear mixed effects (NLME) and quantitative systems pharmacology (QSP). NLME methods like FOCEI and LaplaceI have an estimation phase for the random effects [1] while QSP methods for virtual populations [2] require successive estimation to characterize the uncertainty of the fitting space. However, the robustness of such techniques can be difficult as nonlinear optimization methods can require tuning multiple hyperparameters and are prone to achieving local minima. Here we introduce a technique with a software implementation which allows for robustly identifying all potential solutions to the inference problem by mixing rational interpolation with computer algebra techniques. Using a homotopy technique, we demonstrate the ability to automatically characterize the structural identifiability of the parameters and give the set of all parameters which solve the inverse problem when there are finitely many solutions.
Methods:
The technique follows 6 phases:
- A rational interpolation is fit to the time series data.
- The ODE system is symbolically differentiated sufficiently many times such that substitution can guarantee a polynomial system. This is guaranteed to be finite by computer algebra results derived for structural identifiability analysis [3,4]. This is done by using the internal functionality of the StructuralIdentifiability.jl package.
- The rational interpolation is differentiated to give derivative estimates for derivative values up to the amount required by the symbolically differentiated ODE.
- The derivative values are replaced into the symbolically extended ODE.
- A homotopy technique is used to solve the resulting polynomial nonlinear system and derive all potential solutions (all local minima) to the inverse problem.
- The resulting parameters are assessed by running the ODE solver at those values to calculate loss values at the local minima to classify the result and derive which minima is the global minimum.
The result gives the parameters at each local optima which corresponds to a structurally provable optima (as noted through structural analysis), where the optima are differentiated by the practical identifiability as measured by the cost of the direct ODE solve.
Results:
- This method is tested on a set of 4 ODE models, including a 3-compartment model and a HIV model, all with less than 10 ODEs with a subset of the observables chosen to give cases with global identifiability and local identifiability.
- The tests against other methods like direct single shooting methods with local sensitivity analysis derivatives, a global eSS method from AMIGO2 [5], and the Nelder-Mead derivative-free local optimization method of IQM [6] demonstrate that the new robust parameter estimation method is able to result in an estimation that is approximately an order of magnitude lower in maximum relative error in comparison to the direct shooting, AMIGO2, and IQM methods. This new parameter estimation method is the only one able to demonstrate <100% relative error in the estimated parameters over all ODE models.
- In cases where global structural identifiability is not possible, the method is demonstrated to accurately give estimates for the number of local optima and the relative optimal points.
- The method is shown to be robust in a way that requires no hyperparameter tuning from the user.
Conclusion:
The parameter estimation methods used today in pharmacometrics require repeated simulation and can be prone to many of the issues inherent in fitting with shooting-based techniques. With this method we demonstrate not only the ability to perform the estimation robustly, but also tell the user if there are potentially multiple solutions. However, such a method still has many limitations:
- It is limited to rational ODEs. All compartment models and models with Michaelis-Menton or Hill function activations fall into this category and therefore it’s not a major restriction for the domain though it is a notable restriction of the method.
More research is needed to assess the robustness of the technique to noise. In particular, noise filtering approaches to improve the rational interpolation may be required with very noisy data.
References:
[1] Bonate, Peter L., and Danny R. Howard. Pharmacokinetics in drug development. Arlington, VA: AAPS Press, 2004.
[2] Rieger, Theodore R., et al. "Improving the generation and selection of virtual populations in quantitative systems pharmacology models." Progress in biophysics and molecular biology 139 (2018)
[3] Hongyu Miao, Xiaohua Xia, Alan S. Perelson, and Hulin Wu, On Identifiability of Nonlinear ODE Models and Applications in Viral Dynamics, SIAM Review, 2011.
[4] Alexey Ovchinnikov, Anand Pillay, Gleb Pogudin, and Thomas Scanlon, Computing all identifiable functions for ODE models, preprint, 2020.
[5] Borisov I, Metelkin E (2020) Confidence intervals by constrained optimization—An algorithm and software package for practical identifiability analysis in systems biology. PLoS Comput Biol 16(12).