Efficient Methods for Parameter Estimation of Ordinary and Partial Differential Equation Models of Viral Hepatitis Kinetics

Abstract

Parameter estimation in mathematical models that are based on differential equations is known to be of fundamental importance. For sophisticated models such as age-structured models that simulate biological agents, parameter estimation that addresses all cases of data points available presents a formidable challenge and efficiency considerations need to be employed in order for the method to become practical. In the case of age-structured models of viral hepatitis dynamics under antiviral treatment that deal with partial differential equations, a fully numerical parameter estimation method was developed that does not require an analytical approximation of the solution to the multiscale model equations, avoiding the necessity to derive the long-term approximation for each model. However, the method is considerably slow because of precision problems in estimating derivatives with respect to the parameters near their boundary values, making it almost impractical for general use. In order to overcome this limitation, two steps have been taken that significantly reduce the running time by orders of magnitude and thereby lead to a practical method. First, constrained optimization is used, letting the user add constraints relating to the boundary values of each parameter before the method is executed. Second, optimization is performed by derivative-free methods, eliminating the need to evaluate expensive numerical derivative approximations. The newly efficient methods that were developed as a result of the above approach are described for hepatitis C virus kinetic models during antiviral therapy. Illustrations are provided using a user-friendly simulator that incorporates the efficient methods for both the ordinary and partial differential equation models.

Keywords: constrained optimization; derivative free optimization; differential equations; multiscale models; parameter estimation; viral hepatitis.

Figure A1.

Illustration of the original simplex method. The points x₀, x₁, and x₂ form the initial simplex. (A) The point $\overline{x}$ is the midpoint of the line joining x₀ and x₁, and $\widehat{x}$ is the reflection of x₂ through this line. If $f(\widehat{x})<f\left(x_1\right),x_2$ is replaced by $\widehat{x}$, shifting the location of the simplex. (B) If $f(\widehat{x})\ge f\left(x_1\right),x_2$ is replaced with (1/2)(x₂ + x₀) and x₁ is replaced with (1/2)(x₁ + x₀), reducing the volume of the simplex.

Figure A2.

Illustration of the Nelder–Mead method. The points x₀, x₁ and x₂ form the initial simplex. The vertex x₂ is replaced with a vertex of the form $x_{\text{new }}=\overline{x}+\theta \left(\overline{x}-x_2\right)$. If $f(\widehat{x})<f\left(x_0\right),\theta =2$, if $f\left(x_0\right)\le f(\widehat{x})<f\left(x_{n-1}\right),\theta =1/2$, and if $f\left(x_{n-1}\right)\le f(\widehat{x}),\theta =-1/2$.

Figure A3.

Minimization of $\widehat{f}$ The candidate vertex x* is computed by minimizing $\widehat{f}(x)$ subject to the constraints ĉ₁, ĉ₁ ≥ 0 within the trust region ∥x − x₀∥₂ ≤ ρ. (top) The region of optimization (green) is the intersection of the trust region ∥x − x₀∥₂ ≤ ρ with the half planes defined by the affine constraints ĉ₁ ≥ 0 and ĉ₂ ≥ 0. (bottom left) The function, $\widehat{f}$, to be minimized is represented graphically by the plane (blue) passing through points (x₀, f(x₀)), (x₁, f(x₁)), (x₂, f(x₂)). (bottom right) The vertex x* is defined to be the point within the region of optimization (green) at which $\widehat{f}$ (blue) is minimized.

Figure A4.

Minimization of $\widehat{M}$ Should the constraints ĉ_i(x) ≥ 0 be inconsistent with one another within the trust region ∥x − x₀∥₂ ≤ _ρ, the candidate vertex x* is chosen to minimize $\widehat{M}≔\text{max}\left\{-{\widehat{c}}_i(x):i=1,\dots ,m\right\}$. (top) The constraints ĉ₁(x) ≥ 0 and ĉ₂(x) ≥ 0 are inconsistent within the region ∥x − x₀∥₂ ≤ ρ. (bottom left) Graphs of the affine functions −ĉ₁(x) (blue) and −ĉ₂(x) (green). (bottom right) The vertex x* is defined to be the point within the trust region (black circle) at which $\widehat{M}$ is minimized.

Figure A5.

Illustration of the new vector x^Δ, generated to improve the shape of the simplex. The vertex x₁ is replaced with either $x^{\Delta }=x_0+\gamma \rho v_{\mathcal{l}}$ or $x^{\Delta }=x_0-\gamma \rho v_{\mathcal{l}}$, whichever point results in a smaller value of $\widehat{\mathrm{\Phi }}$.

Figure A6.

Biphasic model fitting example with data taken from [36] of a patient who was treated with mavyret. The LSF method (default) is recommended for use.

Figure A7.

Biphasic model fitting example with data taken from [36] of a patient who was treated with epclusa. In this particular case, COBYLA was selected instead of LSF and succeeded to yield a fit.

Figure A8.

Fitting the parameters c and ρ of the multiscale model to generated data points using the LSF method.

Figure A9.

Fitting the parameters c and ρ of the multiscale model to generated data points using the COBYLA method.

Figure 1.

A flowchart of our constrained damped Gauss–Newton method.

Figure 2.

Start fit that emanates from data of a patient reported in [48]. The fitting curve corresponds to default parameters before fitting with our methods. The multiscale model is used.

Figure 3.

End fit using Gauss–Newton (LSF) that emanates from data of a patient reported in [48].

Figure 4.

End fit using COBLYA that emanates from data of a patient reported in [48].

Figure 5.

Comparison between the line fits of different methods inside the simulator window for the retrieved data points of patient HD that was reported in [48].

Figure 6.

Comparison between the line fits of different methods for the retrieved data points of patient HD that was reported in [48].