J Pharmacokinet Pharmacodyn DOI 10.1007/s10928-017-9527-z

ORIGINAL PAPER

Exploring inductive linearization for pharmacokinetic– pharmacodynamic systems of nonlinear ordinary differential equations Chihiro Hasegawa1,2

•

Stephen B. Duffull1

Received: 11 December 2016 / Accepted: 22 May 2017  Springer Science+Business Media New York 2017

Abstract Pharmacokinetic–pharmacodynamic systems are often expressed with nonlinear ordinary differential equations (ODEs). While there are numerous methods to solve such ODEs these methods generally rely on time-stepping solutions (e.g. Runge–Kutta) which need to be matched to the characteristics of the problem at hand. The primary aim of this study was to explore the performance of an inductive approximation which iteratively converts nonlinear ODEs to linear time-varying systems which can then be solved algebraically or numerically. The inductive approximation is applied to three examples, a simple nonlinear pharmacokinetic model with Michaelis–Menten elimination (E1), an integrated glucose–insulin model and an HIV viral load model with recursive feedback systems (E2 and E3, respectively). The secondary aim of this study was to explore the potential advantages of analytically solving linearized ODEs with two examples, again E3 with stiff differential equations and a turnover model of luteinizing hormone with a surge function (E4). The inductive linearization coupled with a matrix exponential solution provided accurate predictions for all examples with comparable solution time to the matched time-stepping solutions for nonlinear ODEs. The time-stepping solutions however did not perform well for E4, particularly

Electronic supplementary material The online version of this article (doi:10.1007/s10928-017-9527-z) contains supplementary material, which is available to authorized users. & Chihiro Hasegawa [email protected] 1

School of Pharmacy, University of Otago, Dunedin, New Zealand

2

Translational Medicine Center, Ono Pharmaceutical Co., Ltd., Osaka, Japan

when the surge was approximated by a square wave. In circumstances when either a linear ODE is particularly desirable or the uncertainty in matching the integrator to the ODE system is of potential risk, then the inductive approximation method coupled with an analytical integration method would be an appropriate alternative. Keywords Nonlinear ordinary differential equations  Recursive systems  Iterative linearization  Analytical approximations  Matrix exponential solution  Proper lumping

Introduction Pharmacokinetic–pharmacodynamic (PKPD) systems are often expressed using ordinary differential equations (ODEs) as a way of directly translating a biological process into simple input–output equations. In some cases these PKPD ODE systems cannot be expressed with exact closed forms due to inclusion of nonlinearity in which either an input or output function is dependent on the response variable itself (e.g. Michaelis–Menten kinetics [1]). These nonlinear systems are usually solved with time-stepping ODE solvers. Time-stepping solvers are a series of extremely valuable generic routines, for example the Runge–Kutta methods, that can be applied to numerically solve systems of ODEs. The benefit of time-stepping solvers is their general lack of dependence on the model structure, their solution depending only on the information from derivatives (i.e. slope against time) of the response variable within the time step. Many modelling software applications contain inbuilt time-stepping solvers, for example ADVAN6 in NONMEM [2], or provide ready access to these solvers (e.g. in R and MATLAB). The generality of their application, i.e.

123

J Pharmacokinet Pharmacodyn

lack of dependence on the underlying model characteristics, does however come at some limitations including accuracy and computation time under some circumstances. To avoid this problem, there is a need to match the timestepping method to the structure of the problem. A common example is the use of stiff time-stepping solvers (e.g. Gear integrator which is also known to be used in ADVAN8 in NONMEM [2]) in systems where there is a large difference in the rate constants among the differential equations. These solvers, although generically applicable may result in large unexpected inaccuracies in model predictions which do not occur with analytical solutions, due to their ‘‘blinkered’’ time-stepping properties. An inductive approximation has been introduced as a technique for generating approximated solutions to nonlinear systems based on iterative linearization [3]. It uses the intermediate solution obtained from the previous iteration to convert nonlinear ODEs to linear time-varying systems which converges to the system solution to any given level of accuracy. It is noted, however, that even after the inductive approximations, exact closed form solutions may not be always available. However, even in such cases, the system can be analytically as well as numerically solved and therefore provide an approximate solution to an exact problem. In this sense time-stepping ODE solvers provide an approximate solution to an approximate problem. One of the potential advantages of the inductive approximation is the availability of the analytical solution (that may be available via the integrating factor formula [4], Laplace transforms, or matrix exponential solutions [2]) to inform the integration process of the ODEs, and hence some of the potential issues of time-stepping solvers can be avoided. Another advantage lies in generating the linearized system which can be used directly in some settings, e.g. scale reduction (‘‘lumping’’) of complex PKPD systems in which exact solutions are only available for linear ODEs [5] or generating derivatives for construction of the Fisher information matrix [6]. The inductive approximation has been tested only through a single example [3]. The primary aim (Aim 1) of this study was to explore the performance of an inductive approximation to convert nonlinear ODEs to linear timevarying systems. The secondary aim (Aim 2) of this study was to explore the potential advantages of analytical methods over time-stepping solutions. For the primary aim, the inductive approximations are applied to three examples, a simple nonlinear pharmacokinetic (PK) model with Michaelis–Menten elimination (E1), an integrated glucose– insulin (IGI) model [7–9] and an HIV viral load model [10–12] with recursive feedback systems (E2 and E3, respectively). The difference between E2 and E3 is that the latter is also an example of a stiff system. For the secondary aim, two examples are considered, again E3 with

123

stiff differential equations and a turnover model of luteinizing hormone (LH) with a surge function (E4) [13]. This paper is divided into the following sections: theory describing the inductive approximation and thereafter methods to analytically obtain the solution; application to examples (E1–E4) and discussion.

Theory Inductive approximation for systems of nonlinear ordinary differential equations The inductive solution has been described elsewhere [3]. A similar concept including its proof of convergence and an estimate of the error has also been introduced elsewhere as the ‘‘Picard iteration method’’ [14]. Here we briefly cover the main concepts of the inductive approximation. We consider how to solve PKPD systems of nonlinear ODEs. These are commonly defined in the form: dy ¼ fðt; yÞ þ Aðt; yÞy; dt

yðt0 Þ ¼ y0

ð1Þ

where t is time, y(t) = (y1(t), …, ym(t))T is an m 9 1 vector of response variables (i.e. m states in the PKPD system), f(t, y) is an m 9 1 vector, and A(t, y) is an m 9 m nonsingular matrix of which each element includes the corresponding parameters (e.g. rate constants). Generally f corresponds to input functions, and A can be both input and output functions with first-order rate constants. Starting with an initial approximation y[0](t), the solution y[n](t) can be inductively obtained in the form: dy½n ¼ fðt; y½n1 Þ þ Aðt; y½n1 Þy½n ; dt

y½n ðt0 Þ ¼ y0

ð2Þ

Here y[n-1](t) is a response vector obtained from the previous iteration and hence considered to be a known quantity. Since f and A now depend on a known quantity, y[n-1](t), rather than on the current value of y[n](t), then Eq. 2 now represents a linear time-varying system which can be defined in the form: dy ¼ fðtÞ þ AðtÞy; dt

yðt0 Þ ¼ y0

ð3Þ

With the inductive approximation, the response vector can be expressed with the convergent property as:     yðtÞ ¼ y½0 ðtÞ þ y½1 ðtÞ  y½0 ðtÞ þ y½2 ðtÞ  y½1 ðtÞ þ ... ð4Þ and hence, with a sufficient n, the inductive approximation (Eq. 2) will form the true solution of which only the timederivative is known exactly (Eq. 1).

J Pharmacokinet Pharmacodyn

After the inductive approximations, exact closed form solutions to PKPD systems are not always available as the system although linear is time-varying. In this time-varying system, either f or A (or both) (from Eq. 3) depend on t. A typical example of this system is a turnover model with drug concentration–time profiles as a driver of pharmacodynamic markers [15]. However, even in such cases, the system can be analytically as well as numerically solved as the approximation of an exact solution. In this study we apply the following analytical methods for comparison purposes: • • •

Integrating factor formula and solution of integrals Laplace transforms with exact/numerical inversion Matrix exponential solution

Analytical methods to the solution of systems of linear time-varying ordinary differential equations Application of an integrating factor formula Since Eq. 2 is a linear system, we can apply the integrating factor formula [4], see also Duffull et al. [3] for a worked example, to write the solution as: 0 t 1 Z   y½n ðtÞ ¼ exp@ A s; y½n1 ðsÞ dsA y0 to

þ

Zt

0

exp@

Zt

1     A v; y½n1 ðvÞ dvAf s; y½n1 ðsÞ ds

s

t0

ð5Þ In some cases one or both of integrals may be solved exactly. Even if the exact solution cannot be obtained, a closed form analytical approximation can be used to estimate the integrals such as that attained with Gaussian quadrature [16] or simply using a linear trapezoidal method. Laplace transforms Laplace transforms are well documented and used in PKPD and can be used to obtain the solution of Eq. 2 [17]. When A is time-invariant, the solution can be written as: 1 y ðtÞ ¼ 2pi ½n

cþi1 Z

   ðsE  AÞ1 y0 þ F s; y½n1 ðsÞ ds

ci1

ð6Þ where E is an identity matrix, and F(s) is the Laplace transform of the original vector f(t). When the inverse Laplace transform (written in Eq. 6) cannot be performed exactly, it can be approximated

accurately using the rapid Numerical Inverse Laplace Transform (NILT) [18]. An additional modification needs to be considered in the case when A is time-varying since Laplace transforms are no longer straightforward to apply. In this work we explore the use of a linear approximation, i.e. to use time-invariant values of A within the range of a specified step size (i.e. t ± h where h is the step size), such that the degree of timevarying behavior between two consecutive time points provides a solution that remains appropriately accurate. For example, A(t) will be used for Laplace transforms within the range (t-h, t ? h), and for the next range (t ? h, t ? 3h), initial values for responses need to be updated with y[n](t ? h) obtained from the last point of the current range (t-h, t ? h). Note that even though this represents a time-stepping solution the solution is dependent on the structural model in order to find solutions at future time points. Hence we believe this method will provide a way of avoiding the pitfalls of generic ‘‘blinkered’’ time-stepping ODE methods. Matrix exponentials Another way of analytically solving the system written in Eq. 2 is to use the matrix exponential (ME) solution. This method is used in ADVAN5/7 in NONMEM for solving the system before parameter estimation or simulation [2]. As is known, the ME solution generally provides an exact prediction only for linear time-invariant systems which is generally defined in the form: dy ¼ f þ Ay; dt

yðt0 Þ ¼ y0

ð7Þ

In this system, f and A are just time-constant, and the ME solution of Eq. 2 with t0 = 0 can be written as: y½n ðtÞ ¼ expðtAÞ y0 þ ðexpðtAÞ  EÞ A1 f

ð8Þ

Although the decomposition of eigenvalues and eigenvectors may be used to calculate the matrix exponential exp(tA) they do not always provide an accurate solution especially when the system is complex and/or the eigenvector matrix is poorly conditioned or not positive definite. In these circumstances a Pade´ approximation of the matrix exponential which does not rely on the condition of the eigenvector matrix can be used and is a versatile method [19]. Similar to the method of Laplace transforms we handle the time-varying nature of Eq. 2 by creating a linear piecewise approximation for small h. Application In this study four examples are considered:

123

J Pharmacokinet Pharmacodyn

• • • •

a simple nonlinear PK model with Michaelis–Menten elimination (E1; to address Aim 1) a standard integrated glucose–insulin (IGI) model [7–9] (E2; to address Aim 1) an HIV viral load model [10–12] (E3; to address Aim 1 and Aim 2) an empirical surge model [13] (E4; to address Aim 2).

Examples E1 to E3 are nonlinear ODE systems where an inductive linearization is used. Example 4 is a time-varying linear system which provides a further example to assess analytical solutions (although not requiring linearization). To address our aims the solution of a numerical timestepping solver with the original nonlinear ODE system (Eq. 1 rather than Eq. 3) is used as a reference (R1) since there are no closed-form exact solutions for these systems of nonlinear ODEs. Figure 1 shows a flow chart to solve nonlinear ODE systems with the inductive approximation. After the inductive approximation, the nonlinear ODE system is converted to a linear time-varying system, and hence either a time-stepping solver or analytical methods described in the previous section can be used (third box from the top in Fig. 1). All of the methods tested are: • • • • •

Original nonlinear ODE system with time-stepping solver (R1) Linearized ODE system with time-stepping solver (N1) Linearized ODE system with integrating factor formula and trapezoidal method for integration (A1) Linearized ODE system with integrating factor formula and Gauss–Legendre quadrature for integration (A2) Linearized ODE system with exact Laplace transforms (A3)

Fig. 1 Flow chart of the inductive approximation. n number of iterations in inductive approximations

Start with

•

Linearized ODE system with numerical inverse Laplace transforms (NILT) (A4) Linearized ODE system with matrix exponential (ME) solution (A5)

•

The methods A1–A5 correspond to the analytical solutions. The method A4 (NILT) is assessed also in the case when the exact Laplace transforms can be performed. Generally the time-stepping solver ode45 in MATLAB is used for the reference R1 and N1. This is a medium order non-stiff differential equation solver based on the 4th order Runge–Kutta technique and is similar to the method applied in ADVAN6 in NONMEM [2]. This method automatically determines the step-size for time-stepping but can be forced to step at additional time points using the tspan option in which the user specifies additional userdefined time points. MATLAB (MATrix LABoratories R2015b, The Math Works Inc) is used through all examples. SAS (Statistical Analysis Software version 9.3, SAS Institute Inc) is partly used to prepare figures. Simple nonlinear PK model with Michaelis–Menten elimination (E1) Here the inductive approximation is applied to the PK model with first-order input and Michaelis–Menten output, i.e. V

dC Vmax ¼ ka Dexpðka tÞ  C; dt Km þ C

Cðt0 Þ ¼ 0

ð9Þ

where C(t0) is the initial concentration of drug in the central compartment with units of mg/L, D is dose (1 mg,

=1

Select inial approximaon for state(s) causing nonlinearity

Solve ODE for

= Is soluon sufficiently accurate? YES

Stop

123

NO

+1 &

Update with previous soluon

J Pharmacokinet Pharmacodyn

single dose) and the parameters, V (volume of distribution), ka (first-order absorption rate constant), Vmax (maximum velocity of the enzyme), and Km (Michaelis–Menten constant defined as the concentration for which the enzyme is half-maximal rate). The parameter values are given in Table 1. Considering the form of Eq. 1, f ðt; CÞ ¼ f ðtÞ ¼ ka Dexpðka tÞ=V;

Aðt; CÞ ¼ 

Vmax =V Km þ C ð10Þ

and hence, with the inductive approximation starting with the naive initial approximation of C[0](t) = 0 (for all t), we can now generate successive approximations inductively which converge to the solution of Eq. 9, i.e. dC ½n Vmax =V ¼ ka Dexpðka tÞ=V  C ½n ðtÞ; dt Km þ C ½n1 ðtÞ C½n ðt0 Þ ¼ 0

ð11Þ

We consider the solution of Eq. 11 with either the method N1 or A1 to A5. The method A2, Gaussian quadrature, has already been applied in the previous work [3]. The same number of points for abscissas in Gaussian quadrature (=5) was also used in this study. For the trapezoidal method (A1), the increment to approximate the integrals was set to 0.1 h which was also used for the fixed time component used in ode45. For Laplace transforms (A3), the step size h shown in the previous section was set to 0.1 h. For NILT (A4), the step size h was also set to 0.1 h, and the number of points in fast Fourier transformations to 28. For the ME solution (A5) the step size h was set to 1 h, and the decomposition of eigenvalues and eigenvectors was used to calculate matrix exponentials. Note these values of h have been chosen empirically to provide sufficiently accurate solutions. Figure 2 (a, top panel) shows the results of the first six iterations with respect to the reference method (R1) when using method A5. The results using other methods are shown in the Supplementary Material. Figure 3 (a, top panel) shows the results with respect to the solution time

Table 1 Parameter values in the simple nonlinear PK model with Michaelis–Menten elimination Parameter

Description

Value

V

Volume of distribution (L)

1

ka

Absorption rate constant (/h)

1

Vmax

Maximum velocity of the enzyme (mg/h)

0.0734

Km

Michaelis–Menten constant (mg/L)

0.3672

for all methods including the reference R1. All of the methods tested showed that as number of iterations increase the inductive approximations approach the solution provided by the reference R1. A sufficient number of iterations was also similar through all methods tested (n = 6). However, the solution time varied in the order of 105. The method N1 provided the final solution (n = 6) with comparable solution time to the reference R1. Although the solution time of methods A1–A4 was longer compared to the reference R1, the method A5 provided the final solution with comparable or even faster solution time. Integrated glucose–insulin (IGI) model (E2) Here the inductive approximation is applied to the IGI model of which the schematic representation is shown in Fig. 4. This model consists of the following seven states: glucose absorption, glucose transit, central glucose, peripheral glucose, glucose effect compartment, central insulin, and insulin effect compartment. All model equations and parameter values can be found in the original paper published by Jauslin et al. [7]. In this IGI model, insulin concentrations are described by a 1-compartment model with endogenous insulin secretion (ISEC) into the central compartment and firstorder elimination (CLI) as: dI CLI ¼ ISEC ðtÞ   IðtÞ; dt VI

Iðt0 Þ ¼ ISS  VI

  GE ðtÞ IPRG ISEC ðtÞ ¼ ðISS  CLI Þ  ð1 þ Sincr  GT ðtÞÞ  GSS ð12Þ where I(t0) is the initial amount of insulin in the central compartment with units of pmol, and the parameters, ISS (steady-state insulin concentration), VI (central volume of distribution), Sincr (slope coefficient contributing to the incretin effect triggered by the amount of glucose in the transit compartment (GT) following glucose absorption), and GSS (steady-state glucose concentration). GE is the concentration in the glucose effect compartment that accounts for the regulation of insulin secretion with the power coefficient IPRG, and is expressed as a function of the amount of glucose in the central compartment (GC) as F1(GC). Thus, the insulin secretion is clearly affected by the nonlinear function of the glucose concentration via an effect compartment (GE), and consequently the solution I is expressed as a function of GC as: IðtÞ ¼ F2 ðt; GC Þ

ð13Þ

Considering the form of Eqs. 1, 12 may be rewritten as:

123

J Pharmacokinet Pharmacodyn

Fig. 2 Inductive approximations coupled with an ME solution (A5) for nonlinear PK model with Michaelis–Menten elimination (E1, a), integrated glucose-insulin model (E2, b, c for insulin and glucose concentrations, respectively), and HIV viral load model (E3, d, e for

123

viral load and CD4 counts, respectively). Vertical dotted lines in b and c represent timing of meal intake. n number of iterations in inductive approximations, ME matrix exponential (Color figure online)

J Pharmacokinet Pharmacodyn

GLUCOSE

INSULIN

MEAL

GA Insulin secreon

Glucose producon

ka

Incren effect

kGE GT

GP

ka Q

kGE

GE I

GC +

CLG

+ +

CLGI

kIE

IE

CLI

kIE

Fig. 4 Schematic representation of the integrated glucose-insulin model. Full arrows indicate flows, and broken arrows indicate control mechanisms (?, stimulating effect). Q, CLG, CLGI, and ka, kinetic parameters of the glucose; CLI, kinetic parameters of the insulin; kGE and kIE, rate constants for the effect compartments

absorption from the gut is described by a first-order process (rate constant ka) via a transit compartment. The amount of glucose in the central compartment (GC) is expressed in the form: dGC Q ¼ ka  GT ðtÞ þ GPROD þ GP ðtÞ Vp dt ðCLG þ CLGI  IE ðtÞ þ QÞ   GC ðtÞ; VG GC ðt0 Þ ¼ GSS  VG Fig. 3 Solution time in inductive approximations for nonlinear PK model with Michaelis–Menten elimination (E1, a), integrated glucose-insulin model (E2, b), and HIV viral load model (E3, c). An unshaded bar with solid lines represents the averaged solution time per iteration. The whole bar together with both unshaded and shaded areas represent the solution time through all iterations. An unshaded bar with dashed lines represents the solution time of the reference R1 (time-stepping solver with the original nonlinear ODE system, no iterations required). NILT numerical inverse Laplace transform, ME matrix exponential

dI ¼ f ðt; GE Þ þ Aðt; IÞ  IðtÞ ¼ f ðt; F1 ðGC ÞÞ þ A  IðtÞ; dt Iðt0 Þ ¼ ISS  VI f ðt; GE Þ ¼ ðISS  CLI Þ  ð1 þ Sincr  GT ðtÞÞ   GE ðtÞ IPRG  ; Aðt; IÞ GSS CLI ¼ VI

GPROD ¼ GSS  ðCLG þ CLGI  ISS Þ

ð15Þ

where GC(t0) is the initial amount of glucose in the central compartment with units of mmol, GP is the amount of glucose in the peripheral compartment, and the parameters, Q (intercompartmental clearance), VP (peripheral volume of distribution), and VG (central volume of distribution). IE is the concentration in the insulin effect compartment that accounts for the regulation of insulin-dependent glucose clearance CLGI, and is expressed as a function of the amount of insulin in the central compartment (I) as F3(I).Considering the form of Eqs. 1, 15 may be rewritten as: dGC ¼ f ðt; GC Þ þ Aðt; IE Þ  GC ðtÞ dt ¼ f ðtÞ þ Aðt; F3 ðIÞÞ  GC ðtÞ; GC ðt0 Þ ¼ GSS  VG

ð14Þ

Meanwhile, glucose concentrations are described by a 2compartment model with endogenous glucose production (GPROD) into the central compartment. Glucose elimination is expressed with an insulin-independent and an insulindependent clearance (CLG and CLGI, respectively). Glucose

f ðt; GC Þ ¼ ka  GT ðtÞ þ GSS  ðCLG þ CLGI  ISS Þ Q þ GP ðtÞ; Aðt; IE Þ Vp ðCLG þ CLGI  IE ðtÞ þ QÞ ¼ VG

ð16Þ

Substituting Eq. 13 into Eq. 16 gives a single differential equation for GC which indirectly depends on itself as:

123

J Pharmacokinet Pharmacodyn

dGC ¼ f ðtÞ þ Aðt; F3 ðF2 ðGC ÞÞÞ  GC ðtÞ; dt GC ðt0 Þ ¼ GSS  VG

ð17Þ

From Eq. 16 the solution of GC is expressed as a function of I as: GC ðtÞ ¼ F4 ðt; IÞ

ð18Þ

and substituting this into Eq. 14 gives a single differential equation for I which indirectly depends on itself as: dI ¼ f ðt; F1 ðF4 ðIÞÞÞ þ A  IðtÞ; dt

Iðt0 Þ ¼ ISS  VI

ð19Þ

From Eqs. 17–19 the IGI model can obviously be classified into the system of nonlinear ODEs with recursive (feedback) systems. With the inductive approximation starting with the naive initial approximations for the specific parts causing the nonlinearity in ODEs, we can now generate successive approximations inductively which converge to the solution for all seven states. When using G[0] E (t) in Eq. 12 and/or I[0] (t) in Eq. 15, the amount of insulin and E glucose in the central compartment is inductively expressed in the form:   dI ½n ½n ¼ ðISS  CLI Þ  1 þ Sincr  GT ðtÞ dt !IPRG ½n1 GE ðtÞ CLI ½n ð20Þ    I ðtÞ; GSS VI I ½n ðt0 Þ ¼ ISS  VI ½n

dGC Q ½n ½n ¼ ka  GT ðtÞ þ GPROD þ GP ðtÞ Vp dt ½n1 ðCLG þ CLGI  IE ðtÞ þ QÞ ½n  GC ðtÞ;  VG ½n GC ðt0 Þ ¼ GSS  VG

ð21Þ

Here the naive initial approximation for concentrations of insulin and glucose in the effect compartment (G[0] E (t) = GSS and I[0] E (t) = ISS for all t) is used for inductive approximations. Other initial approximations could be [0] chosen including G[0] C (t) in Eq. 14 and/or I (t) in Eq. 16. Differential equations for remaining states are the same with the original equations except the use of y[n] to denote the nth iteration rather than y for responses. We consider the solution of all 7 states in the IGI model with inductive approximations and either N1 or A1–A5 applied in the previous example. For the trapezoidal method (A1), the increment to approximate the integrals was set to 0.01 h, and as an additional naı¨ve initial approximation G[0] P (t) = GSS  VP (for all t) was also used in Eq. 21, to avoid the complexity in approximating the integrals. For Gaussian quadrature (A2), the number of points for abscissas was set to 5. For Laplace transforms

123

(A3), the step size h shown in the previous section was set to 0.01 h. For NILT (A4), the step size h was also set to 0.01 h, and the number of points in fast Fourier transformations to 28. For the ME solution (A5), the step size h was set to 0.01 h, and a Pade´ approximation was used to calculate the matrix exponential as the eigenvector matrix was not positive definite when applying a decomposition to the IGI model. Simulations were conducted based on the original study design where the subject underwent a 24-h meal test study with scheduled meals at 8 a.m., 2 p.m., and 8 p.m. and snacks at 11 a.m., 5 p.m., and 11 p.m. The detail of the study design can be found in the original paper [7]. Figure 2 (b and c, middle panels) shows the results of the first 15 iterations with respect to the reference method (R1) for insulin and glucose concentrations when using method A5. The results using other methods are shown in the Supplementary Material. Figure 3 (b, middle panel) shows the results with respect to the solution time for all methods except for the quadrature method (A2). With the method A2 (Gaussian quadrature) the solution could not be obtained within a month with n = 4, indicating that this method may not be practical for some examples. With all of the methods tested the inductive approximations approached the solution provided by the reference R1 as the number of iterations increased. With the trapezoidal method (A1) the inductive approximation approached the solution with fewer iterations compared to other methods. In this case the method N1 provided the final solution (n = 15) with comparable solution time to the reference R1. The solution times of analytical methods were however longer by 1 or 2 orders of magnitude, albeit the absolute times remained fast (approx. 100 s). Out of all analytical methods tested, the matrix exponential method (A5) provided the final solution with the shortest solution time. HIV viral load model (E3) Here the inductive approximation is applied to the HIV viral load model (latent dynamic model) of which the schematic representation is shown in Fig. 5. This model consists of the following five states: target non-infected CD4 cells (TNI), latently infected CD4 cells (TL), actively infected CD4 cells (TA), infectious viruses (VIN), and noninfectious viruses (VNI). Please note that here ‘‘V’’ represents the viral load with units of copies/mL rather than the volume of distribution. The differential equations in this system are defined in the form: dTNI ¼ k  ð1  INHRTI Þ  c  TNI  VIN  dNI  TNI ; dt dA  dV  ðaL þ dL Þ TNI ðt0 Þ ¼ c  p  ðaL þ fr  dL Þ

ð22Þ

J Pharmacokinet Pharmacodyn

VNI Protease inhibitor

λ

dV

p γ

VIN

dV

fr

TNI

&

dNI

Reverse transcriptase inhibitor

1-fr

TA αL

dA

½n

TL

dL Fig. 5 Schematic representation of the HIV viral load model. Parameter definitions can be found in the text

dTL ¼ ð1  fr Þ  ð1  INHRTI Þ  c  TNI  VIN  aL  TL  dL dt  TL ; ð1  fr Þ  c  TNI ðt0 Þ  VIN ðt0 Þ TL ðt0 Þ ¼ aL þ dL ð23Þ dTA ¼ fr  ð1  INHRTI Þ  c  TNI  VIN þ aL  TL  dA  TA ; dt dV  VIN ðt0 Þ TA ðt0 Þ ¼ p ð24Þ dVIN ¼ ð1  INHPI Þ  p  TA  dV  VIN ; dt k  dNI  TNI ðt0 Þ VIN ðt0 Þ ¼ c  TNI ðt0 Þ

ð25Þ

dVNI ¼ INHPI  p  TA  dV  VNI ; dt

ð26Þ

VNI ðt0 Þ ¼ 0

In addition to the nonlinearity, there is a large difference in the rate constants among the differential equations, e.g. the value between dNI (=0.0085/day) and dV (=30/day) differs in the order of 104. Thus this model is also classified into the system of stiff ODEs. The inductive approximation started with the naive initial approximation only for the infectious viruses as V[0] IN (t) = VIN(t0) (for all t). Successive approximations can now be generated inductively which converge to the solution for all states. As an instance, the target non-infected CD4 cells are inductively expressed in the form:

where k is the activation rate constant of TNI; dNI the death rate constant of TNI; c the infection rate constant of TNI; fr the fraction of TNI which becomes TA; dA the death rate constant of TA; dL the death rate constant of TL; aL the reactivation rate constant of TL; p the viral production rate constant of TA; dV the death rate constant of virus; INHRTI the drug effect by the reverse transcriptase inhibitor (RTI); INHPI the drug effect by the protease inhibitor (PI). INHRTI and INHPI are set to 0.90 and 0.75, respectively. The other parameter values can be found in the original paper published by Lavielle et al. [10]. This HIV viral load model is also classified into the system of nonlinear ODEs as the infection rate of target non-infected CD4 cells (second term in the right side of Eq. 22) and the production rate of the other CD4 cells are expressed with the second-order term of the two states, non-infected CD4 cells (TNI) and infectious viruses (VIN).

dTNI ½n ½n1 ¼ k  ð1  INHRTI Þ  c  TNI ðtÞ  VIN ðtÞ  dNI dt ½n  TNI ðtÞ; dA  dV  ðaL þ dL Þ ½n TNI ðt0 Þ ¼ c  p  ðaL þ fr  dL Þ ð27Þ Differential equations for remaining states are obtained by replacing VIN in the right side in Eqs. 23 and 24 with V[n-1] , and by using y[n] to denote the nth iteration rather IN than y for responses. We consider the solution of all five states in the HIV viral load model with inductive approximations. Here we focus only on the method N1, and one of the analytical methods, A5 (matrix exponentials) based on the accuracy and shorter solution time obtained in the previous two examples. As a numerical time-stepping solver ode15s in MATLAB was used for the reference R1. This is a stiff differential equation solver based on the Gear integrator technique, and is expected to be suitable for solving this ‘‘stiff’’ system. Previous runs with ode45 resulted in similar predictions (see the Supplementary Material) but in run times that were three orders of magnitude longer than ode15s. The stiff solver was also selected for the method N1. For the ME solution (A5), the step size h was set to 10 days, and a Pade´ approximation was used to calculate the matrix exponentials. Simulations were conducted based on the original study design where the subject was followed 1 year after treatment initiation. The detail of the study design can be found in the original paper [10]. Figure 2 (d and e, bottom panels) shows the results of the first 18 iterations with respect to the reference method (R1) for the viral load (VIN ? VNI) and CD4 counts (TNI ? TL ? TA) when using method A5. The results using the method N1 are shown in the Supplementary Material. The final solution for CD4 counts was obtained with fewer iterations compared to the viral load, even when n = 6. Figure 3 (c, bottom panel) shows the results with respect to the solution time for all methods (R1, N1, and A5). The method N1 as well as A5 provided the final solution (n = 18) with comparable solution time to the reference R1.

123

J Pharmacokinet Pharmacodyn

Turnover model of luteinizing hormone (LH) with surge function (E4) Here we explore the potential advantages of analytical methods compared with time-stepping ODE solvers to solve the linear time-varying system. This is equivalent to solving a nonlinear ODE after inductive linearization. Here we consider a turnover model of LH which is a linear timevarying system, i.e. no inductive approximation is required in this system. The key component of this turnover model is the existence of surge function to describe the preovulatory LH surge in women, and the drug (cetrorelix) causes a delay in this surge and consequently ovulation. In the model plasma concentrations of cetrorelix are described by a 2-compartment model with first-order absorption and elimination. In addition, the profile of LH is captured by a turnover model with (i) plasma concentration–time profiles of cetrorelix as a driver of LH suppression, and (ii) concentration–time profiles of cetrorelix in the effect compartment as a driver of delaying the surge in LH. This turnover model is defined in the form:   dLH Cc ¼ Rin  1   ð1 þ surgeÞ  kout  LH; dt IC50 þ Cc Rin LHðt0 Þ ¼ kout 00  14 11 Emax Ce t  T þ 0 EC50 þCe B A þ1C surge ¼ SA  @@ ð28Þ A SW where LH(t0) is the initial concentration of LH with units of mIU/mL, Cc the plasma concentration of cetrorelix, Rin the zero-order release rate constant of LH, kout the first-order elimination rate constant of LH, IC50 the plasma concentration of cetrorelix causing 50% of maximum suppression, Ce the plasma concentration of cetrorelix in the effect compartment, SA the surge amplitude, T0 the LH surge day in the absence of drug, Emax the maximum delay, EC50 the plasma concentration of cetrorelix in the effect compartment causing 50% of the maximum delay, and SW the surge width. The parameter values can be found in the original paper published by Nagaraja et al. [13]. Here the parameter values after the multiple administration of cetrorelix 0.5 mg are used. Based on the accuracy and shorter solution time obtained in the previous examples, we only focus on the method A5 (without inductive approximation). For the ME solution (A5), the step size h was set to 1 h, and a Pade´ approximation was used to calculate matrix exponentials. Simulations were conducted based on the original study design where multiple doses were administered once a day from day 3 to day 16 (14 doses). The detail of the study design can be found in the original paper [13].

123

Figure 6 (a, top panel) shows the results with respect to the accuracy for LH. The reference R1 as well as the method A5 provided an equivalently accurate solution. The solution time was also comparable between the reference R1 and the method A5 (0.48 and 0.10 s, respectively). We also applied the similar but partially reformed following model in which the surge was approximated by a square wave: dLH ¼ Rin  ð1 þ surgeÞ  kout  LH; dt

LHðt0 Þ ¼

surge ¼ SA  HðT0  SW  t  T0 þ SWÞ

Rin kout ð29Þ

where H(.) is an indicator function which takes 1 when the condition shown within brackets meets (0 for the other conditions). In this reformed model the drug effect is not considered (suppose as a ‘‘placebo’’ arm) to make the situation simple. Figure 6 (b, bottom panel) shows the results

Fig. 6 Solutions of LH turnover model with surge function (E4). a original model; b reformed model. LH luteinizing hormone (Color figure online)

J Pharmacokinet Pharmacodyn

with respect to the accuracy in this system. While the method A5 provided the exact solution, the reference R1 could not find the surge (even when the timing of the surge was a required time in the time-stepping solution).

Discussion In this study we explored the performance of an inductive approximation which iteratively transforms systems of nonlinear ODEs to linear time-varying. After the inductive approximation, the system can then be solved algebraically or numerically, and hence we explored the application of possible numerical or analytical methods (N1 or A1–A5) for solving linear time-varying systems. Analytical methods coupled with inductive approximations provided accurate predictions in the examples considered. The setting required for each method (e.g. increment in the trapezoidal method to approximate the integrals in Eq. 5) was chosen so that each method could efficiently but still accurately provide the solution. For example, with the trapezoidal method a smaller increment (0.01, with units of hours) was required for the example E2 compared to the simpler case E1 (0.1, with units of hours). This implies that the setting for each method depends on the complexity of system, especially on how each state changes over time. Out of the analytical methods tested, a matrix exponential (ME) solution (A5) provided the best results from the view point of accuracy and runtimes for all examples. The ME solution also provides the benefit of being directly applicable to future model manipulation, for example the scale reduction (‘‘lumping’’) of complicated PKPD systems. An example illustrating the link between lumping and the ME solution is given in Appendix 1. In contrast, the solution time of exact Laplace transforms (A3) was generally 100 times longer, whereas numerical inversion of the Laplace transforms, NILT (A4) was both accurate and had relatively short solution times making this a feasible option of the analytical methods evaluated. Focusing on accuracy, all of these analytical methods (A3–A5) performed well even with the additional necessity to use time-invariant coefficients in a time-stepping manner (see ‘‘Theory’’ Sect.). Gauss–Legendre quadrature (A2) however could not provide an accurate solution within a practical solution time frame in the example E2. Even in the example E1, the method A2 provided the final solution with the longest solution time out of all the methods tested. This method conditions the scale of the quadrature abscissas on the previous iteration’s response (y[n-1](s) and y[n-1](v) in Eq. 5) which in some cases cannot be obtained directly, rather than on the original time scale (y[n-1](t) in Eq. 5) which can be obtained directly from the previous iteration.

It should also be noted that we need to obtain several y[n-1](s) at each t, according to the number of points for abscissa. This makes the calculation computationally intensive, and hence the solution cannot always be obtained within a realistic timeframe, especially in the complex system with multiple states as with example E2. In contrast to Gaussian quadrature, the trapezoidal method (A1) however directly uses y[n-1](t), and hence can provide an update in the prediction within a shorter solution timeframe. A1 also performed well from the view point of accuracy although the approximation of integrals (in Eq. 5) by the trapezoidal method can be very dependent on the resolution of the step-size. We also saw with A1 that the inductive approximation approached the solution with fewer iterations compared with the other analytical methods for example E2. The whole solution time was however not shorter compared to other methods, such as a matrix exponential solution (A5). Although even with A5 runtimes were in some cases longer compared to numerical time stepping solvers, an environmental change for computation (e.g. parallel computing) may streamline the process to obtain a solution within each (small) step size. We see in this work that the number of iterations required for the examples varied from six (E1 and E3 for CD4 count) to 18 (E3 for viral load). Whereas in the example E1 which is the simplest situation tested n = 6 provided the sufficient solution, n = 15 was required in the example E2. In the stiff system E3, as might be expected, a sufficient number of iterations differed between response variables, i.e. fewer required for CD4 count (n = 6) than the HIV viral load (n = 18). This may be due to the large (direct) drug effect for CD4 count (INHRTI = 90%). We can see from Eq. 27 that the contribution of the second term of the right hand side can be negligible when the value of INHRTI is large enough. In this case the inductive part in the second term (V[n-1] (t)) does not drastically IN affect the profile of CD4 count and hence the required number of iterations for CD4 count is lower than for the HIV viral load. Meanwhile, through the examples, the method to be used after inductive approximations (N1 or A1–A5) mostly did not affect the sufficient number of iterations (see Fig. 2, except the method A1 and A2). This implies that the sufficient number of iterations depends on the complexity of model structure and parameter values, rather than the method to be applied after inductive approximations. The inductive approximations however converge to the reference solution (time-stepping solvers in the original nonlinear ODE systems) for all examples, irrespective of complexity (e.g. E2) or stiffness of the system (e.g. E3). Even when no closed-form exact solutions are available, in most cases a test run of the solution of a numerical time-stepping solver with the original nonlinear ODE system can be used as a reference. An

123

J Pharmacokinet Pharmacodyn

alternative method could involve assessing the accuracy and convergence of the inductive approximations by running another step (i.e. n ? 1) and comparing the differences. Convergence may also be assessed by computing the difference in the vector of responses between successive iterations. Although from the view point of accuracy and solution time the method N1 provided comparable results to the best analytical method A5 (using matrix exponentials), the solution time was dependent on the choice of time-stepping methods based on the stiffness of the problem in the example E3 (ode15s and ode45). This clearly implies that there is a need to match the method to the structure of the problem irrespective of whether time-stepping method is used from the outset or an inductive linearization followed by analytical integration. A clear benefit of linearization tagged with an analytical solution was evident in evaluation of the empirical surge model where time-stepping solvers could not capture the change of response in the square wave example E4. This is a difficult test for a timestepping integrator as no information would be available from the derivatives of the response variable over time and hence the possible change in the gradient not anticipated. This issue is not expected when using analytical methods since these methods are conditioned on the model response. Although the approximation applied may be an extreme case, it can be appropriate in some cases where an external (non-biological) operation exists, e.g. PK model with dialysis for the renally eliminated drug. In conclusion, the inductive linearization coupled with a matrix exponential solution provided accurate predictions from all examples with comparable solution time to the matched time-stepping solutions for nonlinear ODEs. The time-stepping solutions however did not perform well for some cases from the view point of accuracy and solution time. In circumstances when either a linear ODE is particularly desirable or the uncertainty in matching the integrator to the ODE system is of potential risk, then the inductive approximation method would be an appropriate alternative. Compliance with ethical standards Conflict of interest The authors declare that they have no conflict of interest.

of the reduced system, thereby forming groups that retain a clear physical interpretation [20]. The process of proper lumping for linear systems provides both the structure for the reduced states and also the parameter values and initial conditions for these reduced states. However, the parameter values for the reduced states are not available when a system is nonlinear [5]. We can use the result of inductive approximations (Eq. 11, n = 6 as a final inductive result) to overcome the issue. Since this system has two states, the matrix form can be written as: dy ¼ KðtÞy; yðt0 Þ ¼ y0 dt 0 1 dy1 dy B dt C ¼ @ dy A; 2 dt dt 0 1   0 ka Vmax =V A  y1 ; KðtÞy ¼ @ k  a y2 ½6 Km þ y2 =V   D y0 ¼ 0

where y is the vector of original states (depot and central compartments) with units of mg, and K(t) the original parameter matrix which generally corresponds to the firstorder coefficient A(t) in Eq. 3. The parameter definition in this system can be found in the example E1. Please note that since y2 [6] is considered to be a known quantity, Eq. 30 now represents a linear time-varying system. Then the system for the lumped state can also be written in the form: d^ y ^ ¼ kðtÞ  y^; dt

We consider E1 as the example and wish to lump the depot and central compartments. Here we use proper lumping which is a special case of lumping where each of the original states contributes to only one of the pseudo-states

123

y^ðt0 Þ ¼ y^0

ð31Þ

where y^ is the lumped state with units of mg, and k^ the lumped parameter. Please note that since here the number of states in the lumped model is one, y^ as well as k^ are scalars rather than vector and matrix. Now the question is how to get k^ as well as y^0 to construct the lumped model. The lumped state itself can be calculated with a lumping matrix M (here a matrix of dimensions 1 9 2) consisting of 0 and/or 1 s in the form: y^ ¼ My  My ¼ ð 1

Appendix 1: an example of illustrating the link between lumping and matrix exponential solution

ð30Þ

1Þ

y1 y2

 ¼ y1 þ y2

ð32Þ

and conversely, the vector of original states y can be rewritten in the form: y ¼ Mþ y^

ð33Þ

where M? is the Moore–Penrose inverse, sometimes referred to as pseudo inverse of M. This matrix is unique,

J Pharmacokinet Pharmacodyn

and is being considered as a good choice for the inverse of M [21]. In the similar manner to Eq. 32, y^0 can be obtained in the form: y^0 ¼ My0 My0 ¼ ð 1

 1Þ 

D 0

 ¼D

ð34Þ

By using the lumping matrix M, Eqs. 30, 33, and 34, 31 can be rewritten in the form: d^ y dy ¼ M ¼ MKðtÞy ¼ MKðtÞMþ y^; dt dt

y^ðt0 Þ ¼ y^0 ¼ D ð35Þ

By comparing this equation with Eq. 31, now we know ^ k can be calculated as a time-varying but known variable MK(t)M?. Finally an ME solution with a linear piece-wise approximation for small h directly provides a solution without incorporating any additional variables in the form: y^ðtÞ ¼ expðt  MKMþ Þ D

ð36Þ

where K is the time-invariant parameter matrix within the range of a specified step size ‘‘t ± h’’ obtained from K(t) (see ‘‘Theory’’ Sect.). Now we know that a common original parameter matrix is the basis of calculation both in lumping and the ME solution. The decomposition of eigenvalues and eigenvectors or a Pade´ approximation can be used to calculate the matrix exponential in Eq. 36.

References 1. Beal SL (1983) Computation of the explicit solution to the Michaelis-Menten equation. J Pharmacokinet Biopharm 11(6):641–657 2. Beal SL (1994) NONMEM users guides. University of California at San Francisco, San Francisco 3. Duffull SB, Hegarty G (2014) An inductive approximation to the solution of systems of nonlinear ordinary differential equations in pharmacokinetics-pharmacodynamics. J Theor Comput Sci 2:119 4. Adams RA (1999) Calculus: a complete course, 4th edn. Addison Wesley, Boston 5. Gulati A, Isbister GK, Duffull SB (2014) Scale reduction of a systems coagulation model with an application to modeling pharmacokinetic-pharmacodynamic data. CPT Pharmacomet Syst Pharmacol 3:e90

6. Green B, Duffull SB (2003) Prospective evaluation of a D-optimal designed population pharmacokinetic study. J Pharmacokinet Pharmacodyn 30(2):145–161 7. Jauslin PM, Frey N, Karlsson MO (2011) Modeling of 24-hour glucose and insulin profiles of patients with type 2 diabetes. J Clin Pharmacol 51(2):153–164 8. Silber HE, Jauslin PM, Frey N, Gieschke R, Simonsson US, Karlsson MO (2007) An integrated model for glucose and insulin regulation in healthy volunteers and type 2 diabetic patients following intravenous glucose provocations. J Clin Pharmacol 47(9):1159–1171 9. Schneck KB, Zhang X, Bauer R, Karlsson MO, Sinha VP (2013) Assessment of glycemic response to an oral glucokinase activator in a proof of concept study: application of a semi-mechanistic, integrated glucose-insulin-glucagon model. J Pharmacokinet Pharmacodyn 40(1):67–80 10. Lavielle M, Samson A, Karina Fermin A, Mentre´ F (2011) Maximum likelihood estimation of long-term HIV dynamic models and antiviral response. Biometrics 67(1):250–259 11. Jacqmin P, McFadyen L, Wade JR (2010) Basic PK/PD principles of drug effects in circular/proliferative systems for disease modelling. J Pharmacokinet Pharmacodyn 37(2):157–177 12. Ro¨shammar D, Simonsson US, Ekvall H, Flamholc L, Ormaasen V, Vesterbacka J, Wallmark E, Ashton M, Gissle´n M (2011) Non-linear mixed effects modeling of antiretroviral drug response after administration of lopinavir, atazanavir and efavirenz containing regimens to treatment-naı¨ve HIV-1 infected patients. J Pharmacokinet Pharmacodyn 38(6):727–742 13. Nagaraja NV, Pechstein B, Erb K, Klipping C, Hermann R, Locher M, Derendorf H (2003) Pharmacokinetic/pharmacodynamic modeling of luteinizing hormone (LH) suppression and LH surge delay by cetrorelix after single and multiple doses in healthy premenopausal women. J Clin Pharmacol 43(3):243–251 14. Youssef IK, El-Arabawy HA (2007) Picard iteration algorithm combined with Gauss-Seidel technique for initial value problems. Appl Math Comput 190:345–355 15. Dayneka NL, Grag V, Jusko WJ (1993) Comparison of four basic models of indirect pharmacodynamic responses. J Pharmacokinet Biopharm 21:457–478 16. Trefethen L (2008) Is gauss quadrature better than ClenshawCurtis? SIAM Rev 50:67–87 17. Abuhelwa AY, Foster DJ, Upton RN (2015) ADVAN-style analytical solutions for common pharmacokinetic models. J Pharmacol Toxicol Methods 73:42–48 18. Purves RD (1995) Accuracy of numerical inversion of Laplace transforms for pharmacokinetic parameter estimation. J Pharm Sci 84(1):71–74 19. Golub GH, Van Loan CF (1996) Matrix Computations, 3rd edn. The Johns Hopkins Press Ltd., London 20. Dokoumetzidis A, Aarons L (2009) Proper lumping in systems biology models. IET Syst Biol 3(1):40–51 21. Li G, Rabitz H (1990) A general-analysis of approximate lumping in chemical-kinetics. Chem Eng Sci 45:977–1002

123