Home | About | Journals | Submit | Contact Us | Français |

**|**Bioengineering (Basel)**|**v.4(2); 2017 June**|**PMC5590451

Formats

Article sections

Authors

Related links

Bioengineering (Basel). 2017 June; 4(2): 40.

Published online 2017 May 4. doi: 10.3390/bioengineering4020040

PMCID: PMC5590451

Gou-Jen Wang, Academic Editor

DST/NRF Centre of Excellence in the Mathematical and Statistical Sciences, School of Computer Science and Applied Mathematics, University of the Witwatersrand, Johannesburg, Private Bag 3, Wits 2050, South Africa; Email: hg.ude.smia@nuesawulo (O.E.); Email: az.ca.stiw@sbocaJ.noryB (B.J.)

Received 2017 March 4; Accepted 2017 April 28.

Copyright © 2017 by the authors.

Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

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.

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].

After a drug is released from its dosage form, the drug is absorbed into the surrounding tissue, the body, or both. As commented on by Shargel et al. [8], the distribution through and elimination of the drug in the body varies for each patient but can be characterized using mathematical models and statistics. Being able to characterize drug distribution and elimination is an important prerequisite to be able to determine or modify the dosing regimens of individuals and groups of patients. From among the three main types of pharmacokinetic (PK) models—compartment, physiologic and statistical moment approach models—compartmentally-based models are known to be a very simple and useful tool in pharmacokinetics [8]. In essence, a compartment model provides a simple way of grouping all the tissues into one or more compartments where drugs move to and from the central or plasma compartment. Assuming a drug is given by intravenous (I.V.) injection and that the drug dissolves (distributes) rapidly in the body fluids, one may employ a one-compartment PK model that can describe the situation as a tank containing a volume of fluid that is rapidly equilibrated with the drug [8]. The concentration of the drug in the tank after a given dose is governed by two parameters, namely the fluid volume of the tank that will dilute the drug, and the elimination rate of the drug per unit of time—these are both assumed constant for a given drug [8]. Simplistic as this model may be with regards to drug distribution and elimination in the human body, a drug’s PK properties can frequently be described via such a one-compartment open model. Assuming such a model, the drug is both added to and eliminated from a central compartment which represents plasma and highly-perfused tissues that rapidly equilibrate with the drug. When an I.V. dose of drug is given, the drug enters directly into the central compartment while elimination of the drug occurs from the central compartment given that the organs involved in drug elimination, primarily the kidney and liver, are well-perfused tissues.

In a two-compartment model, the drug can move between the central or plasma compartment to and from the peripheral or tissue compartment [8]. Although the peripheral compartment does not represent a specific tissue, the mass balance accounts for the drug present in all the tissues. Knowing the parameters of either the one- or two-compartment model, one can estimate the amount of drug left in the body and the amount of drug eliminated from the body at any time. A drug that follows the pharmacokinetics of a two-compartment model does not equilibrate rapidly throughout the body, as is assumed for a one-compartment model [9]. In the former, the drug distributes into two compartments, the central compartment and the peripheral compartment. The central compartment represents the blood, extracellular fluid, and highly-perfused tissues. The drug distributes rapidly and uniformly in the central compartment. The peripheral compartment contains tissues in which the drug equilibrates more slowly. Drug transfer between the two compartments is assumed to take place by first-order processes. We will be considering this model in our work.

The simplicity and flexibility of the compartment model is the principal reason for its wide application. A major advantage of such models is that the time course of the drug in the body may be monitored quantitatively with a limited amount of data [10]. Furthermore, these models account accurately for the mass balance of the drug in the body and the amount of drug eliminated (Mass balance includes the drug in the plasma, the drug in the tissue pool, and the amount of drug eliminated after dosage administration.). Compartment models have successfully been applied for the prediction of the pharmacokinetics of the drug and the development of dosage regimens. Moreover, compartment models are very useful in relating plasma drug levels to pharmacodynamic (PD) and toxic effects in the body [10]. Underlying physiologic mechanisms can also be obtained via such a model through model testing of the data. Thus, compartment analysis may lead to a more accurate description of the underlying physiologic processes and the kinetics involved. In clinical PK literature, drug data comparisons are based on compartment models and the easy tabulation of important parameters accomplished [10].

In practice, PK models seldom consider all the rate processes ongoing in the body [9]. Due to the complexity models which incorporate such information pose, simplifying assumptions are often made so that solutions may be obtained. Traditional PK models, being simplified mathematical 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.

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:

$$f(x+h)=f(x)+h{f}^{\prime}(x)+\frac{{h}^{2}}{2!}{f}^{\prime \prime}(x)\cdots $$

such that,

$${f}^{\prime}(x)=\frac{f(x+h)-f(x)}{h}+O(h).$$

Under the assumption that $h\to 0$ then,

$${f}^{\prime}(x)\approx \frac{f(x+h)-f(x)}{h}\phantom{\rule{2.em}{0ex}}\mathrm{or}\phantom{\rule{2.em}{0ex}}{f}^{\prime}({x}_{n})\approx \frac{{f}_{n+1}-{f}_{n}}{h}.$$

(1)

Similarly, if we expand $f(x)$ at $f(x-h)$ then,

$$f(x-h)=f(x)-h{f}^{\prime}(x)+\frac{{h}^{2}}{2!}{f}^{\prime \prime}(x)\cdots $$

such that,

$${f}^{\prime}(x)=\frac{f(x)-f(x-h)}{h}+O(h),$$

which, if $h\to 0$, gives,

$${f}^{\prime}(x)=\frac{f(x)-f(x-h)}{h}\phantom{\rule{2.em}{0ex}}\mathrm{or}\phantom{\rule{2.em}{0ex}}{f}^{\prime}({x}_{n})\approx \frac{{f}_{n}-{f}_{n-1}}{h}.$$

(2)

Subtracting the Taylor series expansion of $f(x-h)$ from $f(x+h)$ provides us with:

$$f(x+h)-f(x-h)=2h{f}^{\prime}(x)+\frac{2{h}^{3}}{3!}{f}^{\prime \prime \prime}(x)+\cdots $$

which, under the assumption that $h\to 0$, provides us with a third approximation for the derivative:

$${f}^{\prime}(x)=\frac{f(x+h)-f(x-h)}{2h}\phantom{\rule{2.em}{0ex}}\mathrm{or}\phantom{\rule{2.em}{0ex}}{f}^{\prime}({x}_{n})\approx \frac{{f}_{n+1}-{f}_{n-1}}{2h}.$$

(3)

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.

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 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:

The orders of the discrete representation of the derivative must be equal to the orders of the corresponding derivatives appearing in the differential equations.

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

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.

*Consider a first-order differential equation of the form:*

$$\begin{array}{c}\hfill \frac{du}{dt}=f(t,u,\lambda ).\end{array}$$

(4)

*The conventional denominator*
$\Delta t$
*of the system in (4) is:*

$$\begin{array}{c}\hfill \frac{du}{dt}\to \frac{{u}_{k+1}-{u}_{k}}{\Delta t},\end{array}$$

(5)

*which is replaced by a nonnegative function*
$\varphi (\Delta t)$
*where,*

$$\begin{array}{c}\hfill \varphi (h)=h+O({h}^{2}).\end{array}$$

(6)

*In the above, we have:*

$$\begin{array}{c}\hfill h=\Delta t,\phantom{\rule{1.em}{0ex}}t\to {t}_{k}=hk,\phantom{\rule{1.em}{0ex}}u(t)\to {u}_{k},\end{array}$$

(7)

*where* k *is an integer. The exact discrete representations of Equation (4) are given as:*

$$\begin{array}{c}\hfill \frac{du}{dt}\to \frac{{u}_{k+1}-\phi {u}_{k}}{\varphi (h)},\end{array}$$

(8)

*where,*

$$\phi (h)=1+O(h).$$

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

*The nonlinear terms that occur in*
$f(t,u,\lambda )$
*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:*

$$\begin{array}{c}\hfill {u}^{2}\approx {u}_{k}{u}_{k+1}\phantom{\rule{1.em}{0ex}}and\phantom{\rule{1.em}{0ex}}{u}^{3}\approx \frac{2{u}_{k+1}^{2}{u}_{k}^{2}}{{u}_{k+1}+{u}_{k}},\end{array}$$

(9)

*Erdogan and Ozis’ approximation of the nonlinear term, which was clearly stated in [27], is given as:*

$$\begin{array}{c}\hfill {u}^{2}\approx {u}_{k}{u}_{k+1}\phantom{\rule{1.em}{0ex}}and\phantom{\rule{1.em}{0ex}}{u}^{3}\approx \frac{1}{2}{u}_{k}^{2}(3{u}_{k+1}-{u}_{k}).\end{array}$$

(10)

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.

*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,\lambda )=f({t}_{2},u,\lambda )$
*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.*

The discrete scheme should not introduce extraneous or spurious solutions.

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.

For a system of linear differential equations ${\underset{\_}{x}}^{\prime}=A\underset{\_}{x}$, 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<\infty $. Thus, by evaluating $A\underset{\_}{x}$ 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 $A\underset{\_}{x}=\underset{\_}{0}$. The behavior of these trajectories around equilibrium points provides further insight into the dynamics of the solution, particularly pertaining to varied parameter values.

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:

$$\begin{array}{c}\hfill \frac{d{c}_{i}}{dt}={\displaystyle \sum _{j=0}^{q}}({k}_{ji}{c}_{j}-{k}_{ij}{c}_{i}),\end{array}$$

(11)

where ${c}_{i}(t)$ is the concentration of the drug in compartment *i*, $i=1,2,\dots ,q$ and the rate constant ${k}_{ij}$ governs the rate of the drug movement from the *i*th compartment to the *j*th-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}=\frac{{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 $\frac{D}{{V}_{p}}$ with the unit $\frac{mg/min}{ml}$ [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:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\frac{dc}{dt}\hfill & ={k}_{21}p-{k}_{12}c-{k}_{10}c+I(t),\hfill \\ \frac{dp}{dt}\hfill & ={k}_{12}c-{k}_{21}p.\hfill \end{array}\right.\end{array}$$

(12)

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:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\frac{dc}{dt}\hfill & ={k}_{21}p-{k}_{12}c-{k}_{10}c,\phantom{\rule{1.em}{0ex}}c(0)=1,\hfill \\ \frac{dp}{dt}\hfill & ={k}_{12}c-{k}_{21}p\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}p(0)=0.\hfill \end{array}\right.\end{array}$$

(13)

If $I(t)={I}_{0}$ then (12) becomes a two-compartment I.V. infusion model, given as:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\frac{dc}{dt}\hfill & ={k}_{21}p-{k}_{12}c-{k}_{10}c+{I}_{0},\phantom{\rule{1.em}{0ex}}c(0)=0,\hfill \\ \frac{dp}{dt}\hfill & ={k}_{12}c-{k}_{21}p\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}p(0)=0.\hfill \end{array}\right.\end{array}$$

(14)

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.

We obtain the analytical solution for the blood compartment $c(t)$ and tissue compartment $p(t)$ to Equation (13) via the Laplace transform $\mathcal{L}$. The Laplace transform is an integral transformation where the linear operator $\mathcal{L}\{f(t)\}$ transforms a function $f(t)$ with $t\in {\mathcal{R}}_{\ge 0}$ from the time domain to a function $f(s)$ with $s\in \mathcal{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 $\mathcal{L}\left\{c\right\}$ and $\mathcal{L}\left\{p\right\}$ respectively, we obtain:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}(s+{k}_{10}+{k}_{12}){Q}_{1}-{k}_{21}{Q}_{2}=1\hfill \\ -{k}_{12}{Q}_{1}+(s+{k}_{21}){Q}_{2}=0,\hfill \end{array}\right.\end{array}$$

(15)

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:

$$\begin{array}{ccc}\hfill {Q}_{2}\left(-{s}^{2}-s{k}_{21}-s{k}_{10}-{k}_{10}{k}_{21}-s{k}_{12}-{k}_{12}{k}_{21}+{k}_{12}{k}_{21}\right)& =& -{k}_{12}\hfill \\ \hfill \left({s}^{2}+({k}_{10}+{k}_{12}+{k}_{21})s+{k}_{10}{k}_{12}\right){Q}_{2}& =& {k}_{12}\hfill \\ \hfill {Q}_{2}& =& \frac{{k}_{12}}{{s}^{2}+({k}_{10}+{k}_{12}+{k}_{21})s+{k}_{10}{k}_{12}}\hfill \end{array}$$

providing,

$$\begin{array}{ccc}\hfill {Q}_{2}& =& \frac{{k}_{12}}{(s+{\lambda}_{1})(s+{\lambda}_{2})}\hfill \end{array}$$

(16)

where,

$${\lambda}_{1,2}=\frac{1}{2}\left({k}_{10}+{k}_{12}+{k}_{21}\pm \sqrt{{({k}_{10}+{k}_{12}+{k}_{21})}^{2}-4{k}_{10}{k}_{21}}\right)$$

(17)

such that,

$$\begin{array}{c}\hfill {\lambda}_{1}+{\lambda}_{2}={k}_{12}+{k}_{21}+{k}_{10},\phantom{\rule{1.em}{0ex}}{\lambda}_{1}{\lambda}_{2}={k}_{21}{k}_{10}.\end{array}$$

To obtain the solution for the peripheral compartment we take inverse Laplace transform of Equation (16) which gives:

$$\begin{array}{ccc}{\mathcal{L}}^{-1}\left\{{Q}_{2}\right\}=p(t)\hfill & =\hfill & {\mathcal{L}}^{-1}\left\{\frac{{k}_{12}}{(s+{\lambda}_{1})(s+{\lambda}_{2})}\right\}\hfill \\ & =\hfill & {k}_{12}{\mathcal{L}}^{-1}\left\{\frac{A}{s+{\lambda}_{1}}+\frac{B}{(s+{\lambda}_{2})}\right\},\hfill \end{array}$$

(18)

where,

$$\begin{array}{c}\hfill A=\frac{{k}_{12}}{{\lambda}_{2}-{\lambda}_{1}},\phantom{\rule{28.45274pt}{0ex}}B=\frac{-{k}_{12}}{{\lambda}_{2}-{\lambda}_{1}}\phantom{\rule{4pt}{0ex}}\end{array}$$

such that,

$$\begin{array}{ccc}\hfill p(t)& =& \frac{{k}_{12}}{{\lambda}_{2}-{\lambda}_{1}}\left(exp(-{\lambda}_{1}t)-exp(-{\lambda}_{2}t)\right).\hfill \end{array}$$

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,

$$\begin{array}{ccc}\hfill {k}_{12}c& =& \frac{dp}{dt}+{k}_{21}p(t)\hfill \\ & =& \frac{d}{dt}\left(\frac{{k}_{12}}{{\lambda}_{2}-{\lambda}_{1}}\left(exp(-{\lambda}_{1}t)-exp(-{\lambda}_{2}t)\right)\right)+{k}_{21}\left(\frac{{k}_{12}}{{\lambda}_{2}-{\lambda}_{1}}\left(exp(-{\lambda}_{1}t)-exp(-{\lambda}_{2}t)\right)\right)\hfill \\ & =& \frac{{k}_{12}}{{\lambda}_{2}-{\lambda}_{1}}\left(-{\lambda}_{1}exp(-{\lambda}_{1}t)+{\lambda}_{2}exp(-{\lambda}_{2}t)\right)+\frac{{k}_{12}{k}_{21}}{{\lambda}_{2}-{\lambda}_{1}}\left(exp(-{\lambda}_{1}t)-exp(-{\lambda}_{2}t)\right)\hfill \end{array}$$

such that,

$$\begin{array}{ccc}\hfill c(t)& =& \frac{{\lambda}_{1}-{k}_{21}}{{\lambda}_{1}-{\lambda}_{2}}exp(-{\lambda}_{1}t)+\frac{{k}_{21}-{\lambda}_{2}}{{\lambda}_{1}-{\lambda}_{2}}exp(-{\lambda}_{2}t).\hfill \end{array}$$

The concentration in the central compartment is hence,

$$\begin{array}{c}\hfill c(t)=Eexp(-{\lambda}_{1}t)+(1-E)exp(-{\lambda}_{2}t)\end{array}$$

(19)

and the concentration of the drug in the peripheral compartment is:

$$\begin{array}{c}\hfill p(t)=\frac{{k}_{12}}{{\lambda}_{2}-{\lambda}_{1}}\left(exp(-{\lambda}_{1}t)-exp(-{\lambda}_{2}t)\right),\end{array}$$

(20)

where $E=\frac{{\lambda}_{1}-{k}_{21}}{{\lambda}_{1}-{\lambda}_{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,

$${k}_{21}p-{k}_{12}c-{k}_{10}c+{I}_{0}=0$$

(21)

and,

$${k}_{12}c-{k}_{21}p=0.$$

(22)

Upon substituting (21) into (22) we obtain the particular solution of (14):

$${q}_{p}=\left(\begin{array}{c}c(t)\\ p(t)\end{array}\right)=\frac{{I}_{0}}{{k}_{10}}\left(\begin{array}{c}1\\ \frac{{k}_{12}}{{k}_{21}}\end{array}\right).$$

(23)

To find the complementary solution to Equation (14) at $I(t)=0$) we consider,

$$\left(\begin{array}{c}{c}^{\prime}\\ {p}^{\prime}\end{array}\right)=\left(\begin{array}{cc}-{k}_{10}-{k}_{12}& {k}_{21}\\ {k}_{12}& -{k}_{21}\end{array}\right)\left(\begin{array}{c}c\\ p\end{array}\right)$$

(24)

providing the corresponding matrix,

$$A=\left(\begin{array}{cc}-{k}_{10}-{k}_{12}& {k}_{21}\\ {k}_{12}& -{k}_{21}\end{array}\right).$$

(25)

We can find a fundamental set of solutions for the matrix by finding the eigenvalues and the corresponding eigenvectors. The matrix *A* has eigenvalues $\lambda $ such that:

$${\lambda}^{2}+({k}_{10}+{k}_{12}+{k}_{21})\lambda +{k}_{10}{k}_{21}=0$$

(26)

which, upon solution, provides the following eigenvalues:

$${\lambda}_{1,2}=\frac{1}{2}\left(-{k}_{10}-{k}_{12}-{k}_{21}\pm \sqrt{\left({k}_{10}+{k}_{12}+{k}_{21}\right){}^{2}-4{k}_{10}{k}_{21}}\right).$$

(27)

Suppose $\underset{\_}{v}=\left(\begin{array}{c}{v}_{1}\\ {v}_{2}\end{array}\right)$ are the eigenvectors corresponding to the eigenvalues, then $(A-\lambda I)\underset{\_}{v}=0$. Hence we have that,

$$\left(\begin{array}{cc}-{k}_{10}-{k}_{12}-{\lambda}_{1,2}& {k}_{21}\\ {k}_{12}& -{k}_{21}-{\lambda}_{1,2}\end{array}\right)\left(\begin{array}{c}{v}_{1}\\ {v}_{2}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$$

or,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}(-{k}_{10}-{k}_{12}-{\lambda}_{1,2}){v}_{1}+{k}_{21}{v}_{2}=0\hfill \\ {k}_{12}{v}_{1}+(-{k}_{21}-{\lambda}_{1,2}){v}_{2}=0.\hfill \end{array}\right.\end{array}$$

(28)

From (28), we have that:

$${v}_{1}=\frac{{k}_{21}+{\lambda}_{1,2}}{{k}_{12}}{v}_{2},$$

and,

$${v}_{2}=\frac{{k}_{10}+{k}_{12}+{\lambda}_{1,2}}{{k}_{21}}{v}_{1}.$$

Suppose ${v}_{2}=1,$ then we have that,

$${v}_{1}=\frac{{k}_{21}+{\lambda}_{1,2}}{{k}_{12}},$$

and thus,

$$\underset{\_}{v}=\left(\begin{array}{c}{v}_{1}\\ {v}_{2}\end{array}\right)=\left(\begin{array}{c}\frac{{k}_{21}+{\lambda}_{1,2}}{{k}_{12}}\\ 1\end{array}\right).$$

Let ${\tau}_{1,2}=\frac{{K}_{21}+{\lambda}_{1,2}}{{K}_{12}},$ then we find the corresponding eigenvectors to the eigenvalues to be:

$$\underset{\_}{v}=\left(\begin{array}{c}{v}_{1}\\ {v}_{2}\end{array}\right)=\left(\begin{array}{c}{\tau}_{1,2}\\ 1\end{array}\right).$$

Substituting the value of ${\lambda}_{1,2}$ from (27), we have:

$${\tau}_{1,2}=\frac{-{k}_{10}-{k}_{12}+{k}_{21}\pm \sqrt{\left({k}_{10}+{k}_{12}+{k}_{21}\right){}^{2}-4{k}_{10}{k}_{21}}}{2{k}_{12}}.$$

We can now write our complementary solution as:

$${q}_{c}(t)={A}_{1}{e}^{{\lambda}_{1}t}\left(\begin{array}{c}{\tau}_{1}\\ 1\end{array}\right)+{A}_{2}{e}^{{\lambda}_{2}t}\left(\begin{array}{c}{\tau}_{2}\\ 1\end{array}\right).$$

(29)

The general solution to (14) is then of the form $q(t)={q}_{p}+{q}_{c}$, hence,

$$q(t)=\left(\begin{array}{c}c(t)\\ p(t)\end{array}\right)=\frac{{I}_{0}}{{k}_{10}}\left(\begin{array}{c}1\\ \frac{{k}_{12}}{{k}_{21}}\end{array}\right)+{A}_{1}{e}^{{\lambda}_{1}t}\left(\begin{array}{c}{\tau}_{1}\\ 1\end{array}\right)+{A}_{2}{e}^{{\lambda}_{2}t}\left(\begin{array}{c}{\tau}_{2}\\ 1\end{array}\right),$$

(30)

giving,

$$c(t)=\frac{{I}_{0}}{{k}_{10}}+{A}_{1}{\tau}_{1}{e}^{{\lambda}_{1}t}+{A}_{2}{\tau}_{2}{e}^{{\lambda}_{2}t}$$

(31)

and,

$$p(t)=\frac{{k}_{12}{I}_{0}}{{k}_{10}{k}_{21}}+{A}_{1}{e}^{{\lambda}_{1}t}+{A}_{2}{e}^{{\lambda}_{2}t}.$$

(32)

Given the initial condition $\left(\begin{array}{c}c(0)\\ p(0)\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),$ and employing Equations (31) and (32) we have:

$$\frac{{I}_{0}}{{k}_{10}}+{A}_{1}{\tau}_{1}+{A}_{2}{\tau}_{2}=0,$$

(33)

$$\frac{{k}_{12}{I}_{0}}{{k}_{10}{k}_{21}}+{A}_{1}+{A}_{2}=0.$$

(34)

With some algebraic manipulation we find the arbitrary constants:

$${A}_{1}=\frac{{I}_{0}}{{k}_{10}}\left(\frac{{k}_{12}{\tau}_{2}-{k}_{21}}{{k}_{21}({\tau}_{1}-{\tau}_{2})}\right)$$

(35)

and,

$${A}_{2}=\frac{{I}_{0}}{{k}_{10}}\left(\frac{{k}_{21}-{k}_{12}{\tau}_{1}}{{k}_{21}({\tau}_{1}-{\tau}_{2})}\right).$$

(36)

Our final solution is then given as:

$$q(t)=\left(\begin{array}{c}c(t)\\ p(t)\end{array}\right)=\frac{{I}_{0}}{{k}_{10}}\left(\begin{array}{c}1\\ \frac{{k}_{12}}{{k}_{21}}\end{array}\right)+\frac{{I}_{0}}{{k}_{10}}\left(\frac{{k}_{12}{\tau}_{2}-{k}_{21}}{{k}_{21}({\tau}_{1}-{\tau}_{2})}\right){e}^{{\lambda}_{1}t}\left(\begin{array}{c}{\tau}_{1}\\ 1\end{array}\right)+\frac{{I}_{0}}{{k}_{10}}\left(\frac{{k}_{21}-{k}_{12}{\tau}_{1}}{{k}_{21}({\tau}_{1}-{\tau}_{2})}\right){e}^{{\lambda}_{2}t}\left(\begin{array}{c}{\tau}_{2}\\ 1\end{array}\right).$$

(37)

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 ${\lambda}_{1}<{\lambda}_{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.

Phase portrait representing $c(t)$ and $p(t)$. The arrows illustrates the trajectories of the system.

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.

Phase portrait representing $c(t)$ and $p(t)$. The arrows illustrates the trajectories of the system.

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.

Phase portrait (**a**) shows the steady state concentration is directly proportional to the infusion rate (${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 **...**

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.

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:

$$\begin{array}{ccc}\hfill \frac{{c}_{n+1}-{c}_{n}}{h}& =& {k}_{21}{p}_{n}-{k}_{12}{c}_{n}-{k}_{10}{c}_{n}\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill \frac{{p}_{n+1}-{p}_{n}}{h}& =& {k}_{12}{c}_{n}-{k}_{21}{p}_{n}.\hfill \end{array}$$

(39)

Via the SFD method the system may be written in explicit form as:

$$\begin{array}{ccc}\hfill {c}_{n+1}& =& {c}_{n}+h({k}_{21}{p}_{n}-{k}_{12}{c}_{n}-{k}_{10}{c}_{n})\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill {p}_{n+1}& =& {p}_{n}+h({k}_{12}{c}_{n}-{k}_{21}{p}_{n}).\hfill \end{array}$$

(41)

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:

$$\begin{array}{ccc}\hfill \frac{{c}_{n+1}-\phi {c}_{n}}{\varphi}& =& {k}_{21}{p}_{n}-{k}_{12}{c}_{n}-{k}_{10}{c}_{n}\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill \frac{{p}_{n+1}-\phi {p}_{n}}{\varphi}& =& {k}_{12}{c}_{n}-{k}_{21}{p}_{n}.\hfill \end{array}$$

(43)

where,

$$\varphi =\frac{{e}^{{\lambda}_{1}h}-{e}^{{\lambda}_{2}h}}{{\lambda}_{1}-{\lambda}_{2}}\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}\mathrm{and}\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}\phi =\frac{{\lambda}_{1}{e}^{{\lambda}_{2}h}-{\lambda}_{2}{e}^{{\lambda}_{1}h}}{{\lambda}_{1}-{\lambda}_{2}},$$

and ${\lambda}_{1}$ and ${\lambda}_{2}$ are given as per Equation (17).

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:

$$\begin{array}{ccc}\hfill \frac{{c}_{n+1}-{c}_{n}}{h}& =& {k}_{21}{p}_{n}-{k}_{12}{c}_{n}-{k}_{10}{c}_{n}+{I}_{0}\hfill \end{array}$$

(44)

$$\begin{array}{ccc}\hfill \frac{{p}_{n+1}-{p}_{n}}{h}& =& {k}_{12}{c}_{n}-{k}_{21}{p}_{n}.\hfill \end{array}$$

(45)

This system may be written in explicit form as:

$$\begin{array}{ccc}\hfill {c}_{n+1}& =& {c}_{n}+h({k}_{21}{p}_{n}-{k}_{12}{c}_{n}-{k}_{10}{c}_{n}+{I}_{0})\hfill \end{array}$$

(46)

$$\begin{array}{ccc}\hfill {p}_{n+1}& =& {p}_{n}+h({k}_{12}{c}_{n}-{k}_{21}{p}_{n}).\hfill \end{array}$$

(47)

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:

$${k}_{10}{k}_{21}\ne 0,$$

(48)

gives the general solutions:

$$\begin{array}{c}\hfill c(t)=\frac{1}{{\lambda}_{1}-{\lambda}_{2}}\left(({k}_{10}{k}_{21}-{\lambda}_{2}){c}_{0}+{k}_{21}{p}_{0}+\frac{{I}_{0}({\lambda}_{2}-{k}_{10}{k}_{21}-{k}_{12})}{{k}_{10}}\right){e}^{{\lambda}_{1}(t-{t}_{0})}+\\ \hfill \frac{1}{{\lambda}_{1}-{\lambda}_{2}}\left(({\lambda}_{1}-{k}_{10}{k}_{21}){c}_{0}+{k}_{21}{p}_{0}+\frac{{I}_{0}({k}_{10}{k}_{21}-{\lambda}_{2}+{k}_{12})}{{k}_{10}}\right){e}^{{\lambda}_{2}(t-{t}_{0})}+\frac{{I}_{0}}{{k}_{10}},\end{array}$$

(49)

$$\begin{array}{c}\hfill p(t)=\frac{{\lambda}_{1}+{k}_{10}+{k}_{12}}{{k}_{21}({\lambda}_{1}-{\lambda}_{2})}\left(({k}_{10}{k}_{21}-{\lambda}_{2}){c}_{0}+{k}_{21}{p}_{0}+\frac{{I}_{0}({\lambda}_{2}-{k}_{10}{k}_{21}-{k}_{12})}{{k}_{10}}\right){e}^{{\lambda}_{1}(t-{t}_{0})}+\\ \hfill \frac{{\lambda}_{2}+{k}_{10}+{k}_{12}}{{k}_{21}({\lambda}_{1}-{\lambda}_{2})}\left(({\lambda}_{1}-{k}_{10}{k}_{21}){c}_{0}+{k}_{21}{p}_{0}+\frac{{I}_{0}({k}_{10}{k}_{21}-{\lambda}_{2}+{k}_{12})}{{k}_{10}}\right){e}^{{\lambda}_{2}(t-{t}_{0})}+\frac{{k}_{12}{I}_{0}}{{k}_{10}{k}_{21}},\end{array}$$

(50)

where ${\lambda}_{12}$ is given in Equation (27).

The NSFD scheme of Equation (14) is obtained by making the following transformation in Equations (49) and (50):

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{t}_{0}\to {t}_{k}=hk,\hfill \\ t\to {t}_{k+1}=h(k+1),\hfill \\ {u}_{0}\to {u}_{k},\hfill \\ u(t)\to {u}_{k+1},\hfill \\ {w}_{0}\to {w}_{k},\hfill \\ w(t)\to {w}_{k+1},\hfill \end{array}\right.\end{array}$$

(51)

giving the following results:

$$\begin{array}{c}\hfill {c}_{k+1}=\frac{1}{{\lambda}_{1}-{\lambda}_{2}}\left(({k}_{10}{k}_{21}-{\lambda}_{2}){c}_{k}+{k}_{21}{p}_{k}+\frac{{I}_{0}({\lambda}_{2}-{k}_{10}{k}_{21}-{k}_{12})}{{k}_{10}}\right)\\ \hfill {e}^{{\lambda}_{1}h}+\frac{1}{{\lambda}_{1}-{\lambda}_{2}}\left(({\lambda}_{1}-{k}_{10}{k}_{21}){c}_{k}+{k}_{21}{p}_{k}+\frac{{I}_{0}({k}_{10}{k}_{21}-{\lambda}_{2}+{k}_{12})}{{k}_{10}}\right){e}^{{\lambda}_{2}h}+\frac{{I}_{0}}{{k}_{10}},\end{array}$$

(52)

$$\begin{array}{c}\hfill {p}_{k+1}=\frac{{\lambda}_{1}+{k}_{10}+{k}_{12}}{{k}_{21}({\lambda}_{1}-{\lambda}_{2})}\left(({k}_{10}{k}_{21}-{\lambda}_{2}){c}_{k}+{k}_{21}{p}_{k}+\frac{{I}_{0}({\lambda}_{2}-{k}_{10}{k}_{21}-{k}_{12})}{{k}_{10}}\right){e}^{{\lambda}_{1}h}+\\ \hfill \frac{{\lambda}_{2}+{k}_{10}+{k}_{12}}{{k}_{21}({\lambda}_{1}-{\lambda}_{2})}\left(({\lambda}_{1}-{k}_{10}{k}_{21}){c}_{k}+{k}_{21}{p}_{k}+\frac{{I}_{0}({k}_{10}{k}_{21}-{\lambda}_{2}+{k}_{12})}{{k}_{10}}\right){e}^{{\lambda}_{2}h}+\frac{{k}_{12}{I}_{0}}{{k}_{10}{k}_{21}}.\end{array}$$

(53)

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:

$${E}_{\mathrm{Absolute}}=|{X}_{\mathrm{Analytical}\phantom{\rule{4.pt}{0ex}}\mathrm{solution}}-{X}_{\mathrm{Numerical}\phantom{\rule{4.pt}{0ex}}\mathrm{solution}}|.$$

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.

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 Table 2 and Table 3, given that the NSFD method produces “exact” solutions. More importantly, we see from Figure 6 and Figure 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 concentration of the drug for Case 1 where $h=0.53333$ in the (**a**) central compartment, and (**b**) peripheral compartment.

The concentration of the drug for Case 1 where $h=5.3333$ in the (**a**) central compartment, and (**b**) peripheral compartment for a large time interval.

The concentration of the drug for Case 2 where $h=10$ in the (**a**) central compartment, and (**b**) peripheral compartment. At the steady state, the rate of drug infusion (${I}_{0}$) is equal to the elimination rate (${k}_{10}$).

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.

Once more we find that for each step size employed the NSFD method outperforms those methods it is compared to—see Table 4 and Table 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\in [0,2000]$; as such in a range of $[0,1]$ the value of $h=10$ represents $\Delta t=0.005$, whereas $h=26.667$ is represented by $\Delta 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.

The concentration of the drug for Case 2 where $h=26.667$ in the (**a**) central compartment, and (**b**) peripheral compartment.

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.

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.

O.E. thanks the Deutscher Akademischer Austauschdienst (DAAD) for a scholarship and the University of Witwatersrand, Johannesburg, for a Postgraduate Merit Award. C.H. and B.J. acknowledge support from the National Research Foundation, South Africa under Grant Nos. 91114 and 94005 respectively. Unless otherwise indicated the views expressed in this paper are my views only and not the views of the DST/NRF Centre of Excellence in the Mathematical and Statistical Sciences.

The following abbreviations are used in this manuscript:

PK | pharmacokinetic |

PD | pharmacodynamic |

I.V. | intravenous |

NSFD | nonstandard finite difference |

SFD | standard finite difference |

Author Contributions

Oluwaseun Egbelowo: conducted all mathematical calculations; structured the necessary code; produced all simulations; wrote the article. Charis Harley: conceptualised the work in terms of context and methodology to be employed; contributed to analysis; contributed to the write up of the article. Byron Jacobs: conceptualised the work in terms of context and methodology to be employed; contributed to analysis; contributed to the write up of the article.

Conflicts of Interest

The authors declare no conflict of interest.

1. Buchanan A. Physiologic effects of the inhalation of ether. Lond. Med. Gaz. 1847;39:715–717.

2. Michealis L., Menten M.L. Die Kinetik der Invertinwirking. Biochem. Z. 1913;49:333–369.

3. Widmark E., Tandberg J. Uber die bedingungen f’tir die Akkumulation Indifferenter Narkoliken Theoretische Bereckerunger. Biochem. Z. 1924;147:358–369.

4. Teorell T. Kinetics of distribution of substances administered to the body. I. The extravascular modes of administration. Arch. Int. Pharmacodyn. Ther. 1937;57:205–225.

5. Holford N.H.G., Sheiner L.B. Kinetics of pharmacologic response. Pharmacol. Ther. 1982;16:143–166. doi: 10.1016/0163-7258(82)90051-1. [PubMed] [Cross Ref]

6. Huang S.-M., Abernethy D.R., Wang Y., Zhao P., Zineh I. The utility of modelling and simulation in drug development and regulatory review. J. Pharm. Sci. 2013;102:2912–2923. doi: 10.1002/jps.23570. [PubMed] [Cross Ref]

7. Bonate P.L. Pharmacokinetic-Pharmacodynamic Modelling and Simulation. Springer; New York, NY, USA: 2011. p. 2.

8. Shargel L., Yu A., Wu-Pong S. Applied Biopharmaceutics & Pharmacokinetics. 6th ed. McGraw-Hill; New York, NY, USA: 2012. Introduction to Biopharmaceutics and Pharmacokinetics; pp. 1–17. Chapter 1.

9. Shargel L., Yu A., Wu-Pong S. Applied Biopharmaceutics & Pharmacokinetics. 6th ed. McGraw-Hill; New York, NY, USA: 2012. Multicompartment Models: Intravenous Bolus Administration; pp. 58–79. Chapter 4.

10. Shargel L., Yu A., Wu-Pong S. Applied Biopharmaceutics & Pharmacokinetics. 6th ed. McGraw-Hill; New York, NY, USA: 2012. Physiologic Pharmacokinetic Models, Mean Residence Time, and Statistical Moment Theory; pp. 658–690. Chapter 22.

11. Beňová M., Gombárska D., Dobrcký B. Using Euler’s and Taylor’s expansion method for solution of non-linear differential equation system in pharmacokinetics. Electr. Rev. 2013;89:259–261.

12. Mickens R.E. Dynamic consistency: A fundamental principle for constructing nonstandard finite difference schemes for differential equations. J. Differ. Equ. Appl. 2005;11:645–653. doi: 10.1080/10236190412331334527. [Cross Ref]

13. Mickens R.E. Application of Nonstandard Finite Difference Schemes. 1st ed. World Scientific; Singapore: 2000.

14. Mickens R.E. Nonstandard finite difference schemes for differential equations. J. Differ. Equ. Appl. 2002;8:823–847. doi: 10.1080/1023619021000000807. [Cross Ref]

15. Mickens R.E. A nonstandard finite difference scheme for the diffusionless Burgers equation with logistic reaction. Math. Comput. Simul. 2003;62:117–124. doi: 10.1016/S0378-4754(02)00180-5. [Cross Ref]

16. Mickens R.E. A numerical integration technique for conservative oscillators combining nonstandard finite-difference methods with a Hamilton’s principle. J. Sound Vib. 2005;285:477–482. doi: 10.1016/j.jsv.2004.09.027. [Cross Ref]

17. Liao C., Ding X. Nonstandard finite difference variational integrators for nonlinear Schrödinger equation with variable coefficients. Adv. Differ. Equ. 2013;2013:1–22. doi: 10.1186/1687-1847-2013-12. [Cross Ref]

18. Arenas A.J., González-Parra G., Chen-Charpentier B.M. A nonstandard numerical scheme of predictor-corrector type for epidemic models. Comput. Math. Appl. 2010;59:3740–3749. doi: 10.1016/j.camwa.2010.04.006. [Cross Ref]

19. González-Parra G., Arenas A.J., Chen-Charpentier B.M. Combination of nonstandard schemes and Richardson’s extrapolation to improve the numerical solution of population models. Math. Comput. Model. 2010;52:1030–1036. doi: 10.1016/j.mcm.2010.03.015. [Cross Ref]

20. Jordan P.M. A nonstandard finite difference scheme for nonlinear heat transfer in a thin finite rod. J. Differ. Equ. Appl. 2003;9:1015–1021. doi: 10.1080/1023619031000146922. [Cross Ref]

21. Malek A. Heat Transfer—Mathematical Modelling, Numerical Methods and Information Technology. InTech; Rijeka, Croatia: 2011. Applications of nonstandard finite difference methods to nonlinear heat transfer problems. [Cross Ref]

22. Introduction to Non-Standard Finite-Difference Modelling. [(accessed on 1 January 2016)]; Available online: https://www.crewes.org/ForOurSponsors/ResearchReports/2006/2006-46.pdf.

23. Mickens R.E. Lie methods in mathematical modelling: Difference equation models of differential equation. Math. Comput. Model. 1988;11:528–530. doi: 10.1016/0895-7177(88)90549-3. [Cross Ref]

24. Mickens R.E. Exact solutions to a finite-difference model of a nonlinear reaction-advection equation: Implications for numerical analysis. J. Differ. Equ. Appl. 2003;9:995–1006. doi: 10.1080/1023619031000146959. [Cross Ref]

25. Mickens R.E. Nonstandard Finite Difference Models of Differential Equations. World Scientific; Singapore: 1994.

26. Anguelov R., Lubuma J.M.-S. Nonstandard finite difference method by nonlocal approximation. Math. Comput. Simul. 2003;61:465–475. doi: 10.1016/S0378-4754(02)00106-4. [Cross Ref]

27. Gurski K.F. A simple construction of nonstandard finite-difference schemes for small nonlinear systems applied to SIR models. Comput. Math. Appl. 2013;66:2165–2177. doi: 10.1016/j.camwa.2013.06.034. [Cross Ref]

28. Péchère J.-C., Péchère M.-M., Dugal R. Clinical Pharmacokinetics of Sisomicin: Two-Compartment Model Analysis of Serum Data after I.V. and I.M. Administration. Eur. J. Clin. Pharmacol. 1976;10:251–256. doi: 10.1007/BF00558337. [Cross Ref]

29. Widder D.V. The Laplace Transform. Princeton University Press; Princeton, NJ, USA: 1966.

30. Koch G. Ph.D. Thesis. University of Konstanz; Konstanz, Germany: 2012. Modeling of Pharmacokinetics and Pharmacodynamics with Application to Cancer and Arthritis.

31. Foster D.M. Noncompartmental versus Compartmental Approaches to Pharmacokinetic Analysis. In: Atkinson A.J., Abernethy D.R., Daniels C.E., Dedrick R.L., Markey S.P., editors. Principles of Clinical Pharmacology. 2nd ed. Academic Press; New York, NY, USA: 2014. pp. 89–105.

32. Wang Z., Kim S., Quinney S.K., Zhou J., Li L. Non-compartment model to compartment model pharmacokinetics transformation meta-analysis–a multivariate nonlinear mixed model. BMC Syst. Biol. 2010;4:S8 doi: 10.1186/1752-0509-4-S1-S8. [PMC free article] [PubMed] [Cross Ref]

Articles from Bioengineering are provided here courtesy of **Multidisciplinary Digital Publishing Institute (MDPI)**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |