Nonstandard Finite Difference Method Applied to a Linear Pharmacokinetics Model

We extend the nonstandard finite difference method of solution to the study of pharmacokinetic–pharmacodynamic models. Pharmacokinetic (PK) models are commonly used to predict drug concentrations that drive controlled intravenous (I.V.) transfers (or infusion and oral transfers) while pharmacokinetic and pharmacodynamic (PD) interaction models are used to provide predictions of drug concentrations affecting the response of these clinical drugs. We structure a nonstandard finite difference (NSFD) scheme for the relevant system of equations which models this pharamcokinetic process. We compare the results obtained to standard methods. The scheme is dynamically consistent and reliable in replicating complex dynamic properties of the relevant continuous models for varying step sizes. This study provides assistance in understanding the long-term behavior of the drug in the system, and validation of the efficiency of the nonstandard finite difference scheme as the method of choice.


Introduction
The first attempt into what led to pharmacokinetics was described by Andrew Buchanan in his work Physiological effects of the inhalation of ether [1] in which he pointed out that for short ether inhalations, the speed of recovery of the ether was related to redistribution of ether in the body. Pharmacokinetics is the science of the kinetics of drug absorption, distribution, and elimination (more precisely relating to excretion and the metabolism). The mathematical representation of this work started with Michaelis and Menten [2] who first developed the now well-known Michaelis-Menten equation to describe enzyme kinetics; later this equation was also used to describe the elimination kinetics of drugs. Widmark and Tandberg [3] in 1924 published equations now known as (a) the one-compartment open model with bolus intravenous injection and multiple doses administered at uniform intervals and (b) the one-compartment open model with constant rate intravenous infusion [3]. Though the full concept of pharmacokinetics was introduced by Teorell [4], Holford and Sheiner [5] defined pharmacokinetics as a branch of pharmacology that employs mathematical models to better understand how drugs are absorbed, distributed, metabolized and excreted by the body. It has been well reported recently that the Food and Drug Administration (FDA) and other drug regulatory agencies have been using modeling and simulation to assist in making informed decisions. Furthermore, pharmaceutical companies are expected to justify their dose, their choice of patient population, and their dosing regimen, not just through clinical trials, but also through modeling and simulation [6,7]. expressions, are based on the assumption of a linear relationship between the dose of a drug and its concentration [11]. In a linear model, these rate coefficients, called k, are assumed to be constant. However, such assumptions regarding the linearity of the model do not necessarily describe the actual physical processes as accurately as a non-linear relationship may. In fact, the non-linearities seen in such models are related to drug absorption, distribution, metabolism and excretion, and the pharmacokinetics of drug action. Another example where simplifying assumptions are made pertains to the number of tissue compartments in a perfusion model. Multi-compartment models were developed to explain the observation that, after a rapid I.V. injection, the plasma level-time curve does not decline linearly as a single, first-order rate process [9]. The plasma level-time curve reflects first-order elimination of the drug from the body only after distribution equilibrium, or plasma drug equilibrium with peripheral tissues occurs. Therefore, while the number of tissue compartments in a perfusion model does vary with the drug, the tissues or organs that have no drug penetration are invariably excluded from consideration. Organs such as the brain, the bones, and other parts of the central nervous system are often excluded, as most drugs have little penetration into these organs [8]. To describe each organ separately with a differential equation would make the model very complex and mathematically difficult to solve. Under such circumstances more sophisticated methods of solution need to be employed.
Our research is aimed at the well-known two-compartment model. We have chosen to consider the two-compartment model in this study, given that the one-compartment model assumes immediate distribution of the drug and attainment of equilibrium throughout the body. Realistically however, very few drugs display these characteristics, and hence we turn to a two-compartment model for consideration. While still a simplification of the actual physics of the problem, it is a well-justified model as commented on above. This model is capable of providing data on rates in and out of specific organs, which is of interest. The work conducted here is done with the aim of introducing a numerical method which may be employed in future research for the solution of models which are non-linear and describe multiple compartments. As such, we propose and illustrate the use of a numerical method of solution, namely the nonstandard finite difference method (NSFD), capable of efficiently obtaining solutions which are not only accurate but maintain the underlying dynamics of the system of equations. This choice of method impacts on whether we are able to consider non-compartment models; the NSFD is not amenable to the simulation of non-compartment models as it provides a meta-analysis of the inter-compartment dynamics, whereas non-compartment models are unable to describe these meta-dynamics and instead conduct parameter estimation of the entire system as a whole through the use of experimental data. The advantage of the NSFD method is the ability to predict the concentration-time profile of a drug when there are alterations in the dosing regimen-this would not be possible were one to consider non-compartment analysis. Another advantage of the NSFD method is that it preserves significant properties of the analogous models and consequently gives reliable numerical results even when analytical solutions are not possible. The standard approaches to multi-compartment models assume linear dynamics over the duration of each time step, whereas the NSFD method assumes exponential dynamics. Thus, in the case of a linear model the NSFD method recovers the model dynamics exactly. This paper illustrates the ability of the NSFD method to solve a two-compartment PK model in a stable and robust fashion, with the ability of being extended to non-linear and/or multi-compartment models.

Finite Difference Method
While the implementation of the NSFD method is the focus of this research, we employ the standard finite difference (SFD) method as a means of comparison. We define finite difference methods as numerical methods used for the solution of differential equations by approximating them with difference equations, in which finite differences approximate the derivatives.
We will introduce the concept of the SFD method from Taylor's theorem, where h is termed the step size between the values of the independent variable x. Were we to increase x by h then according to Taylor's theorem we could create a Taylor series expansion: such that, Under the assumption that h → 0 then, such that, Subtracting the Taylor series expansion of f (x − h) from f (x + h) provides us with: which, under the assumption that h → 0, provides us with a third approximation for the derivative: Equations (1)-(3) are known as forward difference, backward difference and central difference approximations to the first order derivative, respectively. The approximation given by Equation (3) has an error of O(h 2 ) and is hence deemed the most accurate of the three approximations provided.

Nonstandard Finite Difference Method
The NSFD schemes developed by Mickens et al. [12][13][14][15][16] were proposed to compensate for the weaknesses of methods such as the SFD methods; numerical instabilities being a prime example. As commented on by Liao and Ding [17], with regard to the positivity, boundedness, and monotonicity of solutions, NSFD schemes have performed better than SFD schemes. Because it is more flexible in its construction, an NSFD scheme can more easily preserve certain properties and structures obeyed by the original equations and can have better dynamical consistency for dynamical problems [17]. The advantages of NSFD methods have been observed when being employed for many numerical applications. González-Parra et al. [18,19] developed some NSFD methods to preserve the positivity condition and population conservation law of biological models. Heat transfer problems have also been considered via this method in Jordan [20] and Malek [21].
The initial foundation of NSFD schemes came from exact finite difference schemes [22]. It is thought that numerical methods that approximate differential systems are expected to be consistent with the original differential systems. Mickens [23] provides a theorem which states that each ordinary differential equation corresponds to an "exact" finite difference scheme. He also introduces the NSFD method for designing schemes that are dynamically consistent with the original differential systems, preserve physical properties, obtain reliable results and require less effort to implement than those obtained via standard methods [24]. In this study, we employ the NSFD method to obtain solutions for the two-compartment I.V. bolus injection and the two-compartment I.V. infusion models. We furthermore, conduct a comparative analysis by also considering an analytical solution, the SFD method and ODE45 in MATLAB (MathWorks, Natick, MA, USA). The purpose of this paper is to obtain an "exact" finite difference scheme for a linear PK model using the procedure of Mickens [14]. In this fashion we hope to prove the degree to which this method may provide meaningful solutions to equations within this context.

NSFD Modeling Fundamental Principles
NSFD methods provide numerical solutions to differential equations by constructing discrete models. They preserve the significant properties of their continuous analogues and consequently give reliable numerical results. The following rules were given by Mickens in [25] for constructing an NSFD scheme: Rule 1. The orders of the discrete representation of the derivative must be equal to the orders of the corresponding derivatives appearing in the differential equations.

Remark 1.
If the order of the discrete representations for derivatives are larger than those occurring in the differential equations, then numerical instabilities will occur.

Rule 2.
Denominator functions for the discrete representations for derivatives must, in general, be expressed in terms of more complicated functions of the step-sizes than those conventionally used.

Remark 2.
Consider a first-order differential equation of the form: The conventional denominator ∆t of the system in (4) is: which is replaced by a nonnegative function φ(∆t) where, In the above, we have: where k is an integer. The exact discrete representations of Equation (4) are given as: where,

Rule 3.
Nonlinear terms must, in general, be modeled by nonlocal discrete representations.

Remark 3.
The nonlinear terms that occur in f (t, u, λ) are approximated in a nonlocal way by a suitable function of several points on the mesh. For example, Mickens' approximation of nonlinear terms [14,26,27] is given as: Erdogan and Ozis' approximation of the nonlinear term, which was clearly stated in [27], is given as: Rule 4. All the special conditions that correspond to either the differential equation and/or its solutions should also correspond to the difference equation and/or its solutions.

Remark 4.
For example, the system described in Equation (8) is called time-invariant if the behavior of the system does not explicitly depend on the absolute time. In other words, if f (t 1 , u, λ) = f (t 2 , u, λ) for any two times t 1 and t 2 , then the system is time-invariant. If the discrete representation of the same model does not also have this property, numerical instabilities may occur.

Rule 5.
The discrete scheme should not introduce extraneous or spurious solutions.

Remark 5.
The discrete representation of the derivative must, in general, converge to the same fixed-point solutions as the corresponding derivative. If it does otherwise, numerical instabilities may occur.

Phase Plane Analysis
For a system of linear differential equations x = Ax, the phase portrait is a representative set of its solutions, plotted as parametric curves (with t as the parameter) on the Cartesian plane tracing the path of each particular solution (x, y) = (x 1 (t), x 2 (t)) where 0 < t < ∞. Thus, by evaluating Ax at a large number of points and plotting the resulting vectors, one obtains a direction field of tangent vectors to solutions of the system of differential equations. This phase portrait is a graphical tool which assists us in visualizing how the solutions of a given system of differential equations would behave as time evolves. In this context, the Cartesian plane where the phase portrait resides is called the phase plane. The parametric curves traced by the solutions are sometimes also called their trajectories. A qualitative understanding of the behavior of the solutions and local stability of the numerical solution near steady states can usually be gained from a direction field. More precise information can be discovered by including in the plot some solution curves or trajectories.
A phase portrait is able to assist us in establishing whether the trajectories of a solution will approach the equilibrium solution as t increases, where the equilibrium solution is obtained when Ax = 0. The behavior of these trajectories around equilibrium points provides further insight into the dynamics of the solution, particularly pertaining to varied parameter values.

Compartmental Models for Pharmacokinetics
PK models deal with the mathematical description of the rates of the drug movement into, within and upon exiting the body. In this work we depict the whole system as two compartments; employing a two-compartment model is often more appropriate mathematically and physiologically when modeling the rate of a drug in the body. These models are used to predict the time course of drugs in the body and to allow maintenance of drug concentration in the therapeutic range by predicting drug levels in each compartment. In this work we will assume that the drug is uniformly distributed within each compartment. The compartments are also considered to be well-stirred and mixing of the drug is assumed to be rapid. Elimination is depicted as occurring in the central compartment, so the drug in the peripheral compartment must transfer back to the central compartment.
For the q-compartment, we have: where c i (t) is the concentration of the drug in compartment i, i = 1, 2, ..., q and the rate constant k ij governs the rate of the drug movement from the ith compartment to the jth-compartment and vice versa. The variables of importance in this problem structure are as follows: • c : Concentration of drug in central compartment.
• p : Concentration of drug in peripheral compartment.
• k 12 : Transfer rate of drug from central to peripheral compartment.
• k 21 : Degradation rate of drug in peripheral compartment.
• k 10 : Clearance rate of drug leaving the central compartment.
The data relating to the pharmacokinetics of sisomicin, a new single component aminoglycoside antibiotic, were obtained from Péchère et al. [28] to test the model in this study. The elimination profile of this antibiotic follows two-compartment model kinetics after I.V. administration. The I.V. bolus injection is structured as the initial condition which has the units of concentration as mg/mL. In turn, the I.V. bolus infusion is administered as a dose, D, with unit mg/min, as shown in Table 1. The equations considered below are focused on obtaining the concentration in each compartment such that C p = A p V p where A p is the amount of drug present in the p th compartment, V p is the volume of the drug in the p th compartment, with C p representing the compartments: c (central) and p (peripheral) [9].
As such, we consider the infusion I to be defined as D V p with the unit mg/min ml [28].  Figure 1 and are obtained from the work by Péchère et al. [28]. Suppose we consider a two-compartment I.V. infusion model. We employ c(t) and p(t) to represent the drug in the central and peripheral compartment, respectively, within the two-compartment model as depicted in Figure 1. The central compartment is identified with the blood while the peripheral compartment describes soft tissue. The absorption phase is omitted because the drug was administered with I.V. infusion at time zero. Then, Equation (11) reads: We will provide both analytical and numerical solutions to this equation for the case: I(t) = 0 and I(t) = I 0 . If I(t) = 0 then (12) becomes a two-compartment I.V. bolus injection model, given as:

Parameter
The values considered are obtained from the work by Péchère et al. [28]. In this work sisomicin kinetic parameters were derived from a two-compartment open model analysis of serum data after I.V. administration. The values for Subjects 1-4 are employed for the dynamical systems analyses; for the numerical solutions obtained we considered the values of Subject 1. We have also considered different values of the dose; realistically these alternate values should correspond to different values of I 0 , however the assumption is that the parameter values are within the correct range to provide us with at least a semblance of the correct dynamics.

Case 1
We obtain the analytical solution for the blood compartment c(t) and tissue compartment p(t) to Equation (13) via the Laplace transform L. The Laplace transform is an integral transformation where the linear operator L{ f (t)} transforms a function f (t) with t ∈ R ≥0 from the time domain to a function f (s) with s ∈ C in an image domain [29,30]. The main advantage of the Laplace transform is that differentiation and integration in the time domain corresponds to simple algebraic operations in the image domain. Transforming Equation (13) and writing Q 1 , Q 2 for L{c} and L{p} respectively, we obtain: upon substitution of the initial conditions and rearranging accordingly. These are simultaneous algebraic equations in Q 1 and Q 2 ; we can solve for these upon employing the method of elimination: s 2 + (k 10 + k 12 + k 21 )s + k 10 k 12 Q 2 = k 12 Q 2 = k 12 s 2 + (k 10 + k 12 + k 21 )s + k 10 k 12 providing, where, such that, To obtain the solution for the peripheral compartment we take inverse Laplace transform of Equation (16) which gives: where, This then allows us to find c(t) by eliminating Q 2 in Equation (15) or by substituting the solution for p(t) in one of the equations given in (13). In this case, it will be simpler to substitute for p(t) in the second equation of (13), giving, such that, The concentration in the central compartment is hence, and the concentration of the drug in the peripheral compartment is: where E = λ 1 −k 21 λ 1 −λ 2 .

Case 2
The general solution for (14) is of the form q(t) = q p + q c where q p is a particular solution and q c is the complementary solution, i.e., the general solution of the homogeneous part of the Equation (14). To obtain a particular solution of Equation (14), we assume c(t) and p(t) to be constant, such that, and, k 12 c − k 21 p = 0.
Upon substituting (21) into (22) we obtain the particular solution of (14): To find the complementary solution to Equation (14) at I(t) = 0 we consider, providing the corresponding matrix, We can find a fundamental set of solutions for the matrix by finding the eigenvalues and the corresponding eigenvectors. The matrix A has eigenvalues λ such that: which, upon solution, provides the following eigenvalues: Suppose v = v 1 v 2 are the eigenvectors corresponding to the eigenvalues, then (A − λI)v = 0.
We can now write our complementary solution as: The general solution to (14) is then of the form q(t) = q p + q c , hence, giving, and, Given the initial condition c(0) p(0) = 0 0 , and employing Equations (31) and (32) we have: With some algebraic manipulation we find the arbitrary constants: and, Our final solution is then given as:

Remarks
As provided in Table 1 the parameter values for k 10 , k 12 and k 21 were obtained from [28]. We also consider the case where k 12 = 0 as a means of comparison; this case indicates that the transfer rate of the drug from the central to the peripheral compartment is zero.
Consider the dynamics presented in Figure 2 for the homogeneous case; we find that the transition from k 12 > 0 to k 12 = 0 does not change the nature of the equilibrium points. In both cases we find that λ 1 < λ 2 < 0, indicating that we have a node which is defined to be asymptotically stable. This has been checked for the range of parameter values given in Table 1, for Subjects 1-4. This indicates that when the transfer rate of the drug from the central to peripheral compartment is positive, or when there is no transfer, the system maintains the same dynamics.  Figure 3 shows that for the transition from k 12 > 0 to k 12 = 0, in the non-homogeneous case, our dynamics once again remain unchanged. We have a nodal sink-defined to be asymptotically stable-at a point where c and p > 0, such that the concentration in both compartments is positive.
In the second instance (k 12 = 0) we obtain a point where p = 0, i.e., only the concentration in the central compartment is not zero. Thus, we further note that our results have shown that the case where k 10 > k 12 = 0 results in a lower steady-state drug concentration.
We also note-as indicated in Figure 4a-that the steady state concentration is directly proportional to the infusion rate I 0 . Hence, a higher value of I 0 will result in a higher steady-state concentration. In Figure 4b we note that a higher elimination rate (k 10 ) will result in a lower steady-state concentration.       (I 0 ). Thus, a higher I 0 will result in a higher steady-state concentration, and (b) shows a higher elimination rate (k 10 ) will result in a lower steady-state concentration.
In terms of an investigation of the impact of I 0 we note that the drug is given by I.V. infusion at a constant zero-order rate to allow accurate control of the drug concentration and to maintain the level of the drug in the body as constant. This allows us to predict the PK action. For physically meaningful results we can only consider constant values of the dose; at different constant values (I 0 ) we obtain a higher infusion rate which results in a higher steady state concentration. We find however, that an increase in I 0 does not have any effect on the dynamics of the system.

Case 1
Before we discuss the implementation of the NSFD scheme we briefly review the SFD schemes of the system of equations given by (13) which are given as: Via the SFD method the system may be written in explicit form as: We now turn to the focus of this section, which is the implementation of the concepts of the NSFD scheme to the linear PK model (13) under discussion [14]. Employing the method discussed in Section 2 we obtain: where, and λ 1 and λ 2 are given as per Equation (17).

Case 2:
As before, we present the SFD schemes of the system of equations given by (14) before discussing the relevant NSFD scheme. These are given as: This system may be written in explicit form as: We now turn to the NSFD implementation as a means of solution. In order to do so we start by considering the system of equations given by (14) which, under the requirement: gives the general solutions: where λ 12 is given in Equation (27). The NSFD scheme of Equation (14) is obtained by making the following transformation in Equations (49) and (50): giving the following results: (53)

Discussion
Due to its simplicity, the compartment model often serves as a "first model" that requires further refinement in order to describe the physiologic and drug distribution processes in the body accurately. While this may be the case, non-compartment models have gained popularity in the literature due to their ability to easily relate experimental data to mathematical models via parameter fitting. Foste [31] points out that upon comparing non-compartment with compartment models, it is not a question of declaring one method better than the other. It is a question of (1) what information is desired from the data; and (2) what is the most appropriate method to obtain this information. Known limitations of non-compartment models include the fact that they do not allow for meta-analysis and the deeper insight that the physiological-based PK models allow for [32]. Another important limitation of non-compartment PK analysis is that it lacks the ability to predict PK profiles when there are alterations in a dosing regimen, which compartment PK methods are capable of [33]. Furthermore, compartmental approaches allow for some level of "physiological" interpretation of what the body does to the drug (i.e., consistency with physiologically reasonable pathways of drug elimination can be maintained). Shargel et al. [8] note that, in comparison to non-compartment models, compartment models are particularly useful when little information is known about the tissues of the compartment(s). In this study non-compartment models are not of interest, not due to their limitations but because non-compartment analysis is based on algebraic equations, whereas compartment models are based on linear or nonlinear differential equations. We instead investigate the use of a method which is designed for the effective solution of differential equations, and as such turn our attention away from non-compartment models.
In this study two-compartment models (the I.V. bolus injection and I.V. infusion models) are considered for simulations. The methods compared are the SFD method, and the in-built function ODE45 in MATLAB, with the focus of the study being the NSFD method. All the simulations were performed by using MATLAB. The absolute error is presented using the formula: The first column, in each table presented, is the number of nodes chosen, the second column is the corresponding step size, the third column is the absolute error between the exact solution and the SFD scheme, the fourth column is the absolute error between the exact solution and ODE45 in MATLAB, and the last column gives the absolute error between the exact solution and the NSFD scheme.
Once again, we remind the reader of the values provided in Table 1 for the parameters: k 10 , k 12 and k 21 [28]. These were derived from a two-compartment open model analysis of serum of sisomicin after I.V. administration.

Case 1: Simulation Results
The results for this case clearly indicate the degree to which the NSFD method outperforms the SFD method, and in fact the in-built function as well-see Figure 5. For each step size employed the absolute error calculated is in fact zero, as per see Tables 2 and 3, given that the NSFD method produces "exact" solutions. More importantly, we see from Figures 6 and 7 that for large step sizes the NSFD method performs exceptionally well, matching the analytical solution for the entire time profile, while the SFD method is not able to do so.
The SFD method deviates from the exact solution the moment the solution changes gradient-see Figure 7. This is a known problem for methods of this nature; we find that as expected the NSFD method does not suffer from this weakness.

Case 2: Simulation Results
Once more we find that for each step size employed the NSFD method outperforms those methods it is compared to-see Tables 4 and 5. In this case in particular, the structuring of the NSFD scheme is shown to be quite complex. The NSFD scheme in this case was obtained via an initial structuring of the exact solution of the problem and making the substitution provided in Equation (51). We note from Figure 8 that the schemes reach the steady state as required. More importantly we observe the oscillations which occur for increasing step sizes h. It must be noted that these step sizes are realistic given that t ∈ [0, 2000]; as such in a range of [0, 1] the value of h = 10 represents ∆t = 0.005, whereas h = 26.667 is represented by ∆t = 0.013. We can easily observe here the degree to which the NSFD method is able to outperform the other two numerical methods employed. The SFD scheme does not match the real dynamics of the system for higher step sizes; we observe oscillations of the SFD method for large h. As the step size increases there is complete blowup of the SFD scheme and the error in ODE45 increases. In terms of an investigation of the impact of I 0 we note that the drug is given by I.V. infusion at a constant zero-order rate to allow accurate control of the drug concentration and to maintain the level of the drug in the body as constant. This allows us predict PK action. For physically meaningful results we can only consider constant values of the dose; at different constant values (I 0 ) we obtain a higher infusion rate which results in a higher steady state concentration. Table 4. The absolute error results for the system of equations given by (14) for c(t) with parameter values k 10 = 0.0094, k 12 = 0.0405, k 21 = 0.0291 and I 0 = 1.  Table 5. The absolute error results for the system of equations given by (14) for p(t) with parameter values k 10 = 0.0094, k 12 = 0.0405, k 21 = 0.0291 and I 0 = 1.

Conclusions
In this article we structured and applied a NSFD numerical scheme to two systems of equations which are both PK models. The first model is an I.V. bolus injection two-compartment model while the second model is an I.V. infusion two-compartment model. The systems are homogeneous and non-homogeneous systems of differential equations, respectively. We present some absolute errors associated with this method upon comparison to the relevant analytical solution, the SFD method and the in-built function ODE45 in MATLAB. The NSFD results indicates superior performance over the SFD method and the in-built MATLAB function, confirming both its stability and robustness.
We conclude that the developed nonstandard schemes preserve the significant properties of their continuous analogues and consequently give reliable numerical results. It was found that the NSFD schemes were stable for large step sizes and display qualitatively accurate results, allowing one to assess and predict the long time behaviour of the drug in the system. Given that the NSFD method produces "exact" numerical schemes, we may conclude that they are well-suited for the solution of pharmacokinetic models.
The flexibility afforded by the NSFD scheme, in terms of its construction, means that the scheme secures consistency with the continuous pharmacokinetic models arising from compartmental or physiological models with respect to the different dynamical characteristics of the systems [17]. As a consequence, we were able to investigate the dynamics of the models and the impact of the parameter value choices employed. We observe asymptotically stable dynamics in each case and note the impact of parameter value choices on the equilibrium states obtained.
It now remains to extend this method to more complex systems of equations, such as non-linear and/or multi-compartment PK models.