Accurate numerical integration of highly stiff pharmacokinetics models using continuous block implicit hybrid one-step collocation methods

Abstract

The differential equations of pharmacokinetic models, obtained from the formulation based on the Fick's perfusion principle and law of conservation of mass action, deals with absorption, distribution, and elimination of drugs by the body systems. In the formulation, the ordinary differential equations obtained may result into highly stiff systems whereby suitable implicit methods have to be employed to solve accurately drug concentration in the body systems, owing to the complex nature of the exact solutions, if they exist. Further, because of the large and increasing interest in the problems of drug kinetics, it is necessary to apply considerably accurate implicit numerical methods in evaluating and in applying them to specific cases. These solutions are to help determine the distribution of drug concentration in different compartments (parts) of the human body network systems as a fundamental step to help in understanding and improving treatments. This is necessary since drug concentration in each compartment is different from one another because of the differences in drug affinity to tissues. From the study, the solution curves obtained show that drug level in the gastrointestinal tract decreases with the passage of time, while drug concentration in the blood increases from zero and reaches its maximum level and then decreases steadily again. We further, compared the model curves with some experimental data plots published in scientific papers. The results obtained from the study can be used extensively for various drug diffusion problems arising in pharmaceutical studies. The presented results widen the applicability of the continuous block implicit hybrid one-step collocation methods to diffusion process which have good number of applications in the drug control, drug dosage and other related problems in pharmaceutical and biomedical industries.

Introduction

There are many routes of drug delivery into the body systems and their corresponding physiochemical characteristics which play very vital roles in therapeutic efficiency (see e.g., Homayun et al. ¹ ). Thus, some researchers gave the criteria for the selection of drugs delivery route, based on the patient's suitability, solubility of the drug, access to the disease site or location and the effectiveness of the drug in dealing with the specific disease. The absorption mechanism and the nature of the drug are the fundamental factors that determine the appropriate systems for achieving the highest bioavailability and effectivity of drugs. It is an established fact that the diffusion of drug from source to the targeted site is an important challenge in the medical sciences. This is because the influence of drug dosage, inflow and outflow in the processing compartment in the human body system may have both favorable and adverse effects in the body. Therefore, drug intake and concentration at the target site through many compartments in biological processes is a big challenge in the research areas, most especially in the medical sciences. Drug taken by patient must be large enough, in order to be effective, however, too high concentration may also have serious side effect.

As a result, there has been a constant need to develop a universal drug delivery method without modification of the drug. The most recent promising positive progress, proposed some solutions of scaling up the microencapsulation technique and parallelization of the microfluidic channels which both significantly increased the systems yield, as well as addressing the other challenges in existence (see Kannan and Przekwas, ² Kannan et al. ³ ). In Ref., ² the authors studied and developed a multiscale absorption and transit (MAT) toolkit which can simulate the dissolution, transport, absorption, distribution metabolism and elimination of orally administered drugs in the human gastrointestinal tract (GIT) at multiple levels. Thus, their MAT is yet to be developed into a commercial product to meet the urgent demands from the pharmacokinetics and biomedical industries. Although, several commercial software packages have been developed by integrating the drug absorption models with physiologically based pharmacokinetics framework and physiology database to provide predictive tools for drug development and pharmacokinetics studies. However, the prediction powers of most of the models or software tools of drug absorption, distribution and elimination are limited in scope and inherently inefficient. Though, through mathematical approximations a simplified set of equations is demonstrated which can identify a bolus dose required to achieve a specified target effect site of drug concentration. ⁴ Observably all the new technological developments are based on mathematical sciences topological structures (MSTS).

Pharmacokinetic models are very good vehicle tools to the optimal frequency and dosing of drug intake. These represent absorption, distribution, decay and excretion of drug. The flow of drugs within the body can be modeled by treating the different parts of the body the drug passes through as compartments. Drug leaves one compartment and enters another at the rate proportional to the concentration of drug present in the first compartment (see Khandy et al. ⁵ ). The rate at which drug moves between compartments can be modeled by the first- order kinetics. The constant of proportionality in this case is mainly determined by the drug, the compartment and the general health of the individual patient. According to some experiments carried out experimentally, the dosage of medicine varies from one individual patient to another depending upon the condition and severity of the disease.^6,7

In this study we considered accurate solutions of two types of models (see Khanday et al., ⁵ Egbelowo, ⁸ Cascone et al., ⁹ and Bailey and Shafer ¹⁰ ), which consist of II and III compartments for drug absorption and circulation through GIT and blood stream. This is to help estimate the drug concentration at different parts of the body within the blood stream.

The drug delivery models combined the mechanism of drug administration in the human body through oral and intravenous routes. For the intravenous drug delivery two different models were considered due to the reversible and irreversible rate constants in the compartments. In the formulation of the pharmacokinetic models, the stiffness phenomenon exhibited in the drug diffusion models needed to be solve using suitable implicit numerical methods with A-stable regions of absolute stability (RAS). This is because mostly in the literature Laplace transform and eigenvalue methods have been applied to calculate for the exact drug concentration in different compartments of blood stream and tissue medium. ⁵

Therefore, the motivation governing the present study is to apply continuous block implicit hybrid one-step collocation methods to the formulated pharmacokinetic models^5,8–10 and compare the solutions obtained in tabular form as well as in graphical form. The convergence analysis elucidates that the continuous block implicit hybrid one-step collocation methods yield exceedingly accurate results, which are validated by comparisons with Laplace transform and eigenvalue methods ⁵ or the nonstandard finite difference method ⁸ or optimized classical fourth-order Runge–Kutta numerical quadrature ¹¹ used in the literature. Ever in literature, to the best of our knowledge, this is the first time we introduce high-order numerical methods for the solution of such problem in the medical sciences. Definition 1.1

Stiff systems:

The linear system

$$y^{\mathrm{\prime }}=Ay+\mathit{\Phi }(x)$$

is said to be stiff if

• Re(λ_t) < 0, t  =  1,…,p

• (ii) $\underset{t=1,\dots ,p}{max}|\mathrm{Re}({\lambda }_t)|><!--\; >\; -->><!--\; >\; -->\underset{t=1,\dots ,p}{min}|\mathrm{Re}({\lambda }_t)|$

where λ_t are the eigenvalues of the matrix A.

The model problems

Pharmacokinetic model which concerns the study of drug dosage taken into the body systems, either through oral or intravenous route, deals with the inflow, concentration and outflow of the medication which are essential features in the medical sciences for efficient treatment of diseases. In the pharmacokinetic model the formulation of differential equations (ordinary or partial) were based on the Fick's perfusion principle and law of conservation of mass action, simply called the balanced law. ⁵ In general, compartment models describe the transfer of a substance in and out of the compartment, as shown in Figures 1 to 4. This can simply be put as, the rate of change of the amount of the substance (drug) in the compartment is equal to the difference between the rate in which the substance enters the compartment and the rate at which it leaves. For example, if we let $y(t)$ denotes the concentration of drug in the compartment of the body at the time t (in hour) then the rate of change of $y(t)$ is defined by

$${\displaystyle \frac{dy}{dt}=\mathrm{Transfer}\mathrm{in}\text{–}\mathrm{Transfer}\mathrm{out}.}$$

(1)

Figure 1.

Physical representation of oral drug administration through the stomach and blood stream.

Figure 2.

Typical process of drug administration through blood and tissue with initial intake.

Figure 3.

Method of drug administration through arterial blood, tissue and venous with initial intake.

Figure 4.

Schematic representation of a three-compartment drug administration model with initial intake.

Based on the above equation we consider four models for drug distribution taken through either oral or intravenous, depending upon the condition of the patient and the severity of the disease. The exact drug concentration in the different compartments of blood stream and tissue medium were usually calculated by the Laplace transform and eigenvalue methods in the literature. ⁵

Model I

Among the various delivery routes, oral delivery has been recognized as the most widely, attractive and commonly employed route of drug administration worldwide (see/, Ensign et al., ¹² Kannan and Przekwas ¹³ ) because of its potential for solid formulation with long shelf life, sustained delivery and intensified immune response. Further, due to its high possibility, advantages and improved patient compliance, convenience, bioavailability (Sastry et al. ¹⁴ ), ease of administration (Savjani et al. ¹⁵ ) and pain avoidance (Dey and Maiti ¹⁶ ), oral drug delivery is the most widely used and the most readily accepted form of drug administration. When drug is orally taken, it dissolves in the stomach (GI), and releases the medication into the GIT, which diffuse into the blood network. From there, the blood stream takes the medications to the targeted site where it has therapeutic effect. Although, oral bioavailability dependent on the drug compound as well as the physiological and anatomical states of the user, its disadvantages, in some cases considered to be less efficient due to so many reasons which include, stomach sensitivity, liver dysfunction, delayed reaction, chemical and enzymatic stability of drugs, etc.

Generally speaking, the problem of absorption and elimination by an organism through its natural processes are, in reality, quite complicated and could be analyzed more easily after (an explicatory) a clear formulation of the general compartment theory is presented. For example, as in (1), let $y(t)$ be a given concentration of drug in the compartment at a given time t (hour), which can be divided into two parts as follows:

$$y(t)=y_1(t)+y_2(t)$$

(2)

where $y_1(t)$ and $y_2(t)$ represent the concentration of drug for the two phases. Let E be the total fraction of the given dose absorbed, such that the amount of the drug absorbed can be expressed as

$$ED=E_{GI}D+E_B{D}_B$$

where D is the given dose, $E_{GI}$ , $E_B$ denote respectively the fraction of dose absorbed in the two phases ( $y_1(t)$ and $y_2(t)$ ) and $D_B$ is the dose present in the second compartment (see Cherruault and Sarin). ¹⁷

The system of equations describing the rate of change in oral drug administration in the body is given pictorially by Figure 1 ⁵ :

$$\begin{array}{c}{\displaystyle \frac{d{y}_1(t)}{dt}=-c_1{y}_1(t),}{y}_1(0)=D,\\ {\displaystyle \frac{d{y}_2(t)}{dt}=c_1{y}_1(t)-c_e{y}_2(t),}{y}_2(0)=0.\end{array}$$

(3)

In the given system of equations $y_1(t)$ and $y_2(t)$ represent the concentration of drug in the compartment1 (GI) and compartment 2, respectively, at the instant of time t (hour). The values $c_1$ and $c_e$ are the exchange parameters involved in the given system. Specifically $c_1$ is the rate of drug in GI tract, while $c_e$ is the clearance constant as shown in Table 1, so that the given system in (3) are linear in nature. These equations represent a single dosage of drug that flows through the stomach tract to the tissue.

Table 1.

Parameter values for the different models of drug concentration.

Model	Parameter	Interpretation	Value
Model I	$D$	Initial concentration of drug dosage	500 units
	$c_1$	Rate constant of drug in GI tract	0.9776 h−1
	$c_e$	Clearance constant	0.2213 h−1
Model II	$D$	Initial concentration of drug dosage	500 units
	$c_b$	Rate constant of drug in blood stream	0.9776 h−1
	$c_t$	Rate constant of drug in tissue	0.3293 h−1
	$c_e$	Clearance constant	0.2213 h−1
Model III	$D$	Initial concentration of drug dosage	500 units
	$c_{ab}$	Rate constant of drug in arterial blood	0.9776 h−1
	$c_t$	Rate constant of drug in tissue	0.3293 h−1
	$c_e$	Clearance constant	0.2213 h−1

Model II

Medication can also be injected directly into the veins, which is referred to as the intravenous administration. The intravenous drug delivery on the other hand, where the medication is injected directly into the veins, the drug immediately enters the circulatory system without encountering any of the physical, chemical or biological barriers, is regarded as the best suited route for emergency purposes (cases); since the absorption of the drug is guaranteed immediately and uniquely provides precise control of the drug dosage and speed of administration for drugs requiring a stringent dosage. However, its disadvantages entail risks of injection, infection by the needle at the site of injection, circulatory overload, phlebitis and thrombosis. For the intravenous drug administration, two different mathematical models are required on the basis of the reversible rate constant and irreversible rate constant as shown in Figure 2. The main exchange compartments are only two in number, the blood stream and the tissue medium. These are the blood stream into which the drug is injected and the tissue medium where the drug has therapeutic effect. In Figure 2, we illustrate the governing system of two ordinary differential equations describing the rate of change of drug concentration with time, obtained from the two compartments. ⁵

Here we denote the concentration of drug in the blood compartment and tissue compartment by $y_b$ and $y_t$ , respectively, so that the system of equations describing the process are given by

$$\begin{array}{cc}{\displaystyle \frac{d{y}_b(t)}{dt}=-(c_b+c_e)y_b(t)+c_t{y}_t(t),}& y_b(0)=D,\\ {\displaystyle \frac{d{y}_t(t)}{dt}=c_b{y}_b(t)-c_t{y}_t(t),}& y_t(0)=0\end{array}$$

(4\rm a)

where D appearing in the first equation denotes the initial concentration of drug injected through the intravenous to the body. The exact solutions of the equations in (4a) are given by

$$\begin{array}{rl}{y}_b(t)& ={\displaystyle \frac{D}{{\lambda }_2-{\lambda }_1}\{(-{\lambda }_1+c_t)\exp (-{\lambda }_{1}t)-(-{\lambda }_2+c_t)\exp (-{\lambda }_{2}t)\}}\\ y_t(t)& ={\displaystyle \frac{D{c}_b}{{\lambda }_1-{\lambda }_2}\{\exp (-{\lambda }_{1}t)-\exp (-{\lambda }_{2}t)\}}\end{array}$$

(4b)

where ${\lambda }_i$ are the real eigenvalues and are calculated as

$${\lambda }_i={\displaystyle \frac{1}{2}\left[-(c_b+c_t+c_e)\pm \sqrt{{(c_b+c_t+c_e)}^2-4c_e{c}_t}\right],\mathrm{where}i=1,2.}$$

(4c)

Model III

Drug administration into the body through intravenous infusion takes the pattern of Figure 3, because blood flow through the cardiovascular system is in unidirectional. Hence, in the third compartment the effective infinite volume of the drug restricted the transport of substance (drug) essentially one way. In fact, the compartment models which describe the transfer of a substance in and out of the compartment, is as shown in Figure 3.

The model presented in Figure 3 is a pictorial representation of the following system of ordinary differential equation, with parameter values given in Table 2.

$$\begin{array}{cc}{\displaystyle \frac{d{y}_{ab}(t)}{dt}=-c_b{y}_{ab}(t),}& y_{ab}(0)=D,\\ {\displaystyle \frac{d{y}_t(t)}{dt}=c_b{y}_{ab}(t)-c_t{y}_t(t),}& y_t(0)=0,\\ {\displaystyle \frac{d{y}_{vb}(t)}{dt}=c_t{y}_t(t)-c_e{y}_{vb}(t),}& y_{vb}(0)=0.\end{array}$$

(5)

Table 2.

Parameter values of the remifentanil pharmacokinetic models.^8,9

Parameter	Interpretation	Value
$c_{12}$	Transfer rate of drug from central to tissue compartment	0.3730 min−1
$c_{21}$	Degradation rate of drug in tissue compartment	0.1030 min−1
$c_{13}$	Transfer rate of drug from central to deep tissue compartment	0.0367 min−1
$c_{31}$	Degradation rate of drug in deep tissue compartment	0.0124 min−1
$c_{10}$	Clearance rate of drug leaving the central compartment	0.1720min−1
$D$	Initial concentration of drug dosage	5.0000 $\mu g.k{g}^{-1}$

Model IV

Figure 4 shows the direct bolus injection into the central compartment^10,11,18 where the exchanges take place with the other two compartments. This will now allow rapid drug delivery and ensures that the entire dose enters the general circulation, which provides immediate drug effect, especially good for emergency purposes (cases).

The transfer of drug into and out of the compartment are governed by the hybrid rate constants, $c_{12}$ , $c_{21}$ , $c_{13}$ and $c_{31}$ . Elimination from the first or central compartment is given by $c_{10}$ and the initial dose of drug (D) is administered into the first or central compartment at time zero (t  =  0). The differential equations governing this drug delivery administration is given as follows:

$$\begin{array}{cc}{\displaystyle \frac{d{y}_1(t)}{dt}=-(c_{10}+c_{12}+c_{13})y_1(t)+c_{21}{y}_2(t)+c_{31}{y}_3(t),}& y_1(0)=D>0,\\ {\displaystyle \frac{d{y}_2(t)}{dt}=c_{12}{y}_1(t)-c_{21}{y}_2(t),}& y_2(0)=0,\\ {\displaystyle \frac{d{y}_3(t)}{dt}=c_{13}{y}_1(t)-c_{31}{y}_3(t),}& y_3(0)=0.\end{array}$$

(6a)

The exact solutions of the system in (6a) are obtained as

$$y_1(t)=D\left({\displaystyle \frac{(-c_{21}{\xi }_1-c_{31}{\xi }_1+{\xi }_1^2+c_{21}{c}_{31})\exp (-{\xi }_{1}t)}{({\xi }_1-{\xi }_2)({\xi }_1-{\xi }_3)}}\right)+D\left({\displaystyle \frac{(c_{21}{\xi }_3+c_{31}{\xi }_3-{\xi }_3^2-c_{21}{c}_{31})\exp (-{\xi }_{3}t)}{({\xi }_1-{\xi }_3)({\xi }_3-{\xi }_2)}}\right)+D\left({\displaystyle \frac{(c_{21}{\xi }_2+c_{31}{\xi }_2-{\xi }_2^2-c_{21}{c}_{31})\exp (-{\xi }_{2}t)}{({\xi }_1-{\xi }_2)({\xi }_2-{\xi }_3)}}\right)$$

(6b)

$$y_2(t)=D{c}_{12}{\displaystyle \frac{(c_{31}-{\xi }_1)\exp (-{\xi }_{1}t)}{({\xi }_1-{\xi }_2)({\xi }_1-{\xi }_3)}-D{c}_{12}{\displaystyle \frac{(c_{31}-{\xi }_2)\exp (-{\xi }_{2}t)}{({\xi }_1-{\xi }_2)({\xi }_2-{\xi }_3)}-D{c}_{12}{\displaystyle \frac{(c_{31}-{\xi }_3)\exp (-{\xi }_{3}t)}{({\xi }_1-{\xi }_3)({\xi }_3-{\xi }_2)}}}}$$

(6c)

$$y_3(t)=D{c}_{13}{\displaystyle \frac{(c_{21}-{\xi }_1)\exp (-{\xi }_{1}t)}{({\xi }_1-{\xi }_2)({\xi }_1-{\xi }_3)}-D{c}_{13}{\displaystyle \frac{(c_{21}-{\xi }_2)\exp (-{\xi }_{2}t)}{({\xi }_1-{\xi }_2)({\xi }_2-{\xi }_3)}-D{c}_{13}{\displaystyle \frac{(c_{21}-{\xi }_3)\exp (-{\xi }_{3}t)}{({\xi }_1-{\xi }_3)({\xi }_3-{\xi }_2)}}}}$$

(6d)

where ${\xi }_1$ , ${\xi }_2$ and ${\xi }_3$ are referred to as eigenvalues of the system (Upton ¹⁸ ).

Numerical computational techniques

Explicit linear multistep or Runge–Kutta methods tend to have undesirable time-step restrictions when applied to convection–diffusion problems. Implicit linear multistep time-discretization schemes for partial differential equations have proved to be very useful in many real and practical application fields, such as engineering, medicine, science and technology. However, they tend to suffer the disadvantages of very poor stability property as the step number increases with accuracy and requiring additional starting values with constant step size from other one-step methods, like Taylor series, Picard method, Runge–Kutta methods, etc. Runge–Kutta methods, which are the best known one-step methods due to their high accuracy and the ability to incorporate some function evaluations at some off-grid points within the interval of interest, are widely used in solving real life problems. However, they require many function evaluations to achieve better accuracy. Hence, block method which combined both the advantages of fewer function evaluations and dense outputs of the implicit linear multistep methods as well as retaining the A-stability and high accuracy of the Runge–Kutta methods are better option for calculating drug concentration in the compartments. The block methods as they are, can be seen as a set of linear hybrid multistep methods which can be applied simultaneously to pharmacokinetics models and then combined to yield approximations with better accuracy, better efficiency and better stability in parallel computing techniques.

In the past, some kind of continuous Runge–Kutta and linear multistep methods have been derived. The continuous methods were derived to provide dense output of ordinary differential equations (see Yakubu et al. ¹⁹ ) where approximate solutions defined not only on the mesh points but also on some off-grid points as well.

In this research paper, we have proposed new efficient techniques as alternative methods to some well-known methods^5,8,11 to solve pharmacokinetic models. The new techniques are called continuous block implicit hybrid one-step collocation methods, developed based on Legendre polynomial nodes that have better stability regions over a wide range of parameters (see Yakubu et al. ²⁰ and Yakubu and Kwami ²¹ ). The motivation for these methods, particularly the use of Gauss family of methods, is that collocation at Gauss points lead to methods which are symmetric and algebraically stable (see e.g. Hairer and Wanner ²² and Burrage and Butcher ²³ ). It was also testified by Ascher and Bader, ²⁴ that the only symmetric algebraically stable collocation methods are those based on Gauss points. Hence, the numerical computation of our study is based on the continuous block implicit hybrid one-step collocation methods obtained from the Gaussian points, with the general form as follows

$$y(x)=\sum _{j=0}^{t-1}{\varphi }_j(x)y_{n+j}+h\sum _{j=0}^{s-1}{\psi }_j(x)f({\overline{x}}_j,y({\overline{x}}_j))$$

(7a)

where ${\varphi }_j(x)$ and ${\psi }_j(x)$ are the continuous coefficients of the method and are also assumed polynomials of the form,

$${\varphi }_j(x)=\sum _{i=0}^{t+m-1}{\varphi }_{j,i+1}{x}^i\mathrm{and}h{\psi }_j(x)=\sum _{i=0}^{s+m-1}h{\psi }_{j,i+1}{x}^i$$

(7b)

where t denotes the number of interpolation points $x_{n+j,}j=0,1,2,\dots ,t-1$ and s denotes the distinct collocation points ${\overline{x}}_{n+j},j=0,1,2,\dots ,s-1$ chosen from the interval $[x_n,x_{n+1}]$ and the step size h in (7) can be a constant or varied in the numerical integration process. Here, y and f are smooth real N-dimensional vector functions. The numerical constants ${\varphi }_{j,i+1}$ and ${\psi }_{j,i+1}$ in (7b) are to be determined. They are selected so that accurate approximation of well-behaved problem can be obtained efficiently. Theorem 3.1

Collocation at the Gaussian points together with the two end points lead to order one less than the superconvergence scheme, while collocation at only one additional end point to the Gaussian points resulted into order two less than the upgraded scheme.

Proof

see the Ref. ²⁵

Case I: From the proposed collocation scheme (7) after some transformation using combination of LU-decomposition and Gauss-elimination algorithms, we obtain the continuous scheme as,

$$y(x)={\varphi }_0(x)y_n+[{\psi }_1(x)f_{n+u}+{\psi }_2(x)f_{n+v}]$$

(8)

where

$$\xi =(x-x_n),{\varphi }_0(x)=1,{\psi }_1(x)=\left[{\displaystyle \frac{{\xi }^2+2vh\xi }{2(v-u)h}}\right],{\psi }_2(x)=\left[{\displaystyle \frac{{\xi }^2-2uh\xi }{2(v-u)h}}\right],$$

and u, v are the zeros of the Legendre polynomial of degree m (see Yakubu et al. ¹⁹ ). To be specific $m=2$ , that is, Legendre polynomial of degree $m=2$ .

Evaluating the continuous scheme (8) at the points $x=x_{n+1}$ , $x=x_{n+u}$ and $x=x_{n+v}$ we obtain the continuous block implicit hybrid one-step collocation method of order 4 with only 2 stages (Gauss method) written in the form:

$$y_{n+1}=y_n+{\displaystyle \frac{h}{2}[f_{n+u}+f_{n+v}],\mathrm{order}p=4,C_5=2.3148\times {10}^{-4}}$$

(9a)

$$y_{n+u}=y_n+{\displaystyle \frac{h}{12}\left[3f_{n+u}+\left(3-2\sqrt{3}\right)f_{n+v}\right],\mathrm{order}p=2,C_3=8.0187\times {10}^{-3}}$$

(9b)

$$y_{n+v}=y_n+{\displaystyle \frac{h}{12}\left[\left(3+2\sqrt{3}\right)f_{n+u}+3f_{n+v}\right],\mathrm{order}p=2,C_3=-8.0187\times {10}^{-3}}$$

(9c)

From the continuous block implicit hybrid one-step collocation method (9), the well-known Butcher's implicit Runge–Kutta method ²⁶ of order four can be recovered with two stages, if desired. Using Maple 17 version, gives the stability polynomial of the continuous block implicit hybrid one-step method as

$$P(z)={\displaystyle \frac{12+6z+z^2}{12-6z+z^2}}$$

The region of absolute stability of the continuous block method is plotted using Mathmatica version 8.0 package and displayed in Figure 5(a).

Figure 5.

Region of absolute stability of the continuous implicit block hybrid one- step collocation methods.

Case II: For m  =  3, the Legendre polynomial clearly have three zeros (roots), then the continuous scheme obtained is given by

$$y(x)={\varphi }_0(x)y_n+[{\psi }_1(x)f_{n+u}+{\psi }_2(x)f_{n+w}+{\psi }_3(x)f_{n+v}]$$

(10)

where

$$\begin{array}{rl}& {\varphi }_0(x)=1,{\psi }_1(x)=\left[{\displaystyle \frac{2{\xi }^3-3(v+u)h{\xi }^2+6vw{h}^2\xi }{6(v-u)(w-u)h^2}}\right]\\ & {\psi }_2(x)=\left[{\displaystyle \frac{-2{\xi }^3+3(u+w)h{\xi }^2-6uw{h}^2\xi }{6(v-u)(w-v)h^2}}\right],\\ & {\psi }_3(x)=\left[{\displaystyle \frac{2{\xi }^3-3(u+v)h{\xi }^2+6uv{h}^2\xi }{6(w-u)(w-v)h^2}}\right]\end{array}$$

and also u, w and v are the zeros of the Legendre polynomial of degree m  =  3.

Evaluating the continuous scheme (10) at the points $x=x_{n+1}$ , $x=x_{n+u}$ , $x=x_{n+w}$ and $x=x_{n+v}$ we obtain the continuous implicit block hybrid one-step collocation method of order 6 with only 3 stages (Gauss method) written in the form:

$$y_{n+1}=y_n+{\displaystyle \frac{h}{18}[5f_{n+u}+8f_{n+w}+5f_{n+v}],\mathrm{order}p=6,C_7=4.9603\times {10}^{-7}}$$

(11a)

$$\begin{array}{rl}{y}_{n+u}& =y_n\\ & +{\displaystyle \frac{h}{180}\left[25f_{n+u}+\left(40-12\sqrt{15}\right)f_{n+w}+\left(25-6\sqrt{15}\right)f_{n+v}\right],}\\ & \mathrm{order}p=3,C_4=-4.1666\times {10}^{-4}\end{array}$$

(11b)

$$\begin{array}{rl}{y}_{n+w}& =y_n\\ & +{\displaystyle \frac{h}{72}\left[\left(10+3\sqrt{15}\right)f_{n+u}+16f_{n+w}+\left(10-3\sqrt{15}\right)f_{n+v}\right],}\\ & \mathrm{order}p=3,C_4=5.2083\times {10}^{-4}\end{array}$$

(11c)

$$\begin{array}{rl}{y}_{n+v}& =y_n\\ & +{\displaystyle \frac{h}{180}\left[\left(25+6\sqrt{15}\right)f_{n+u}+\left(40+12\sqrt{15}\right)f_{n+w}+25f_{n+v}\right],}\\ & \mathrm{order}p=3,C_4=4.1666\times {10}^{-4}\end{array}$$

(11d)

From the continuous block implicit hybrid one-step collocation method, the well-known Butcher's implicit Runge–Kutta method ²⁶ of order 6 with 3 stages can be recovered, if needed. From the Maple software package we obtain the stability polynomial of the continuous block method as

$$p(z)={\displaystyle \frac{7290+3645z+729z^2+37z^3}{120-60z+12z^2-z^3}}$$

The region of absolute stability of the block method is shown in Figure 5(b) Theorem 3.2

The Gauss method which uses three collocation points only, at $x_{n+u}$ , $x_{n+w}$ and $x_{n+v}$ has order m  =  3 at those points and order 2 m  =  6 at $x_{n+1}$ .

The Radau method which uses four collocation points only, at $x_{n+u}$ , $x_{n+w}$ , $x_{n+v}$ and $x_{n+1}$ has order m  +  1  =  5 at the Gaussian points $x_{n+u}$ , $x_{n+w}$ , $x_{n+v}$ and order m  =  6 at $x_{n+1}$ .

The Lobatto method which uses five collocation points only, at $x_n$ , $x_{n+u}$ , $x_{n+w}$ , $x_{n+v}$ and $x_{n+1}$ has order m  +  2  =  5 at $x_{n+u}$ , $x_{n+w}$ , $x_{n+v}$ and order 2 m  =  6 at $x_{n+1}$ .

Proof

see Ref. ²⁵

Results and discussion

In this section, we consider the solutions of the formulated pharmacokinetic models for the estimation of drug distribution through both oral and intravenous cases. We solve the models using continuous implicit block hybrid one-step collocation methods of orders 4 and 6 with only 2 and 3 stages (Gauss methods) presented in (9) and (11), respectively. The simulations were performed using written codes in MATLAB software, which help to obtain many function evaluations to estimate accurately drug concentration in each of the compartment. We provide four different numerical examples (models) to illustrate the theoretical estimates in nature, but the results established lead to realistic outcome if compared and verified empirically which give the new numerical schemes a better role (position) among numerical methods for the pharmacokinetics models. The numerical results obtained from the computations are displayed graphically as well as in tabular form. The Gauss-points collocation methods produce very fascinating, efficient and consistent results which depict exactly the physical interpretation of the systems of drug delivery. As shown in both the Table of values, as well as the graphical representations, there are better agreement between the computed results and the exact solutions. This further confirmed the efficiency and accuracy of the Gauss collocation methods as stated in the literature by some eminent authors.^22–24

In example 1 or model I, we first treat the whole body system as a single compartment and the graphs obtained of the simulations are shown in Figure 6. In the figure, the drug concentration behavior versus time in the single compartment at different value of the rate constants is shown with an initial dosage of drug as 500 units. It is observed from the graphs of the figure that the drug concentration gradually decreases with respect to time. As clearly depicted from the efficiency curves of the figure, the more the rate constants, the faster the drug absorption within the body system, while the value of the rate constant is small the absorption is less with some portion of the drug retained within the body system which may cause side effects. Practically, this process is shown in the graphs of the figure when the rate constant is 0.9776 hour the drug absorption is very fast but when the rate constant is taken as 0.2213 hour the drug absorption is very slow retaining some portion of the drug within the body system as residue which may cause side effects. Hence to avoid such process two compartments models are considered as shown in Figure 2, with the elimination rate constants playing very important role in the drug absorption. As shown in the plots of Figure 7 the behavior of drug concentration in both the stomach (GI) and blood stream, with the initial dosage of 500 units of drug in the stomach, decreases with the passage of time. As indicated by the plots of Figure 8 the drug concentration in the blood stream increases from zero level to the maximum height and then start to decrease steadily again.

Figure 6.

Drug distribution pattern within the body system in a single compartment of administration with different rate constants (model I). (a) Graph of Method(9). (b) Graph of Method (11).

Figure 7.

Drug concentration pattern within the body through blood and tissue with different rate constants. (a) Graph of the Ode 45 Method. (b) Graph of Method (11).

Figure 8.

Process of drug concentration pattern within the body through arterial blood, tissue and venous with different rate constants. (a) Graph of Method(9). (b) Graph of Method (11).

In the tabular form we compared the results obtained using the new derived continuous implicit block hybrid one-step collocation methods with the results from the Ode45 of MATLAB on which package do exist. From the Table of values, we obtain quite satisfactory results by applying the continuous block implicit hybrid one-step collocation methods (9) and (11). It is evident from the results of Table 3(a) and (b) and Table 4 that with the increase in the order of the approximation method, the accuracy also increases as shown in the performance of methods (9) and (11) with orders 4 and 6, respectively. It is also clear from the results of the tables that the absolute errors are indications that the derived method (11) is generally very accurate for solving the stiff pharmacokinetic models compared with results from the Ode45 code of MATLAB. As clearly shown the absolute errors of the derived method in (11) are very small compared to the errors produced from the results of the Ode codes of MATLAB, which shows that our methods have better position (role) in the numerical computational methods.

Table 3a.

Absolute errors of numerical solutions of model I or Example 1.

$x_i$	$y_i$	Method (9) $|y(x)-y_n(x)|$	Method (11) $|y(x)-y_n(x)|$
5	$y_1$	1.815857331166626 × 10−08	5.684341886080802 × 10−14
	$y_2$	4.766434358316474 × 10−09	1.136868377216160 × 10−13
	$y_3$	8.378719940083101 × 10−11	0
	$y_4$	1.159605744760484 × 10−11	0
50	$y_1$	7.722655936959200 × 10−08	3.694822225952521 × 10−13
	$y_2$	2.607896476547467 × 10−08	8.526512829121202 × 10−14
	$y_3$	7.188987183326390 × 10−10	5.684341886080802 × 10−14
	$y_4$	1.118110048992094 × 10−10	1.136868377216160 × 10−13
250	$y_1$	3.562606876172936 × 10−09	1.465494392505207 × 10−14
	$y_2$	3.686047911344303 × 10−09	1.332267629550188 × 10−14
	$y_3$	7.497362730646273 × 10−10	1.421085471520200 × 10−14
	$y_4$	1.961950601980789 × 10−10	0
500	$y_1$	2.000628829268436 × 10−11	7.372574772901430 × 10−17
	$y_2$	8.390688144288561 × 10−11	3.191891195797325 × 10−16
	$y_3$	2.074767024851099 × 10−10	5.329070518200751 × 10−15
	$y_4$	1.039524022417027 × 10−10	3.552713678800501 × 10−14

Table 3b.

Absolute errors of numerical solutions of model I or Example 1.

$x_i$	$y_i$	Ode45Method $|y(x)-y_n(x)|$	Method (11) $|y(x)-y_n(x)|$
5	$y_1$	7.389644451905042 × 10−13	0
	$y_2$	4.224176564093796 × 10−12	5.126565838509123 × 10−12
50	$y_1$	5.115907697472721 × 10−12	1.136868377216160 × 10−13
	$y_2$	3.905142875737511 × 10−11	4.556000021693762 × 10−11
150	$y_1$	4.845901457883883 × 10−12	5.684341886080802 × 10−14
	$y_2$	6.593836587853730 × 10−11	7.202061169664376 × 10−11
250	$y_1$	2.462030579408747 × 10−12	7.105427357601002 × 10−15
	$y_2$	6.355094228638336 × 10−11	6.650680006714538 × 10−11
350	$y_1$	1.060485033121950 × 10−12	7.105427357601002 × 10−15
	$y_2$	5.317701834428590 × 10−11	5.451283868751489 × 10−11
500	$y_1$	2.613464999967619 × 10−13	1.776356839400251 × 10−15
	$y_2$	3.740296961041167 × 10−11	3.763034328585491 × 10−11

Table 4.

Absolute errors of numerical solutions of model III or Example 3.

$x$	$y_i$	Ode45Method $|y(x)-y_n(x)|$	Method (11) $|y(x)-y_n(x)|$
5	$y_1$	7.389644451905042 × 10−13	0
	$y_2$	5.126565838509123 × 10−12	6.171063660076470 × 10−12
	$y_3$	4.184430579812215 × 10−13	3.558264793923627 × 10−14
50	$y_1$	5.115907697472721 × 10−12	1.136868377216160 × 10−13
	$y_2$	4.527578312263358 × 10−11	5.289280124998186 × 10−11
	$y_3$	3.019806626980426 × 10−12	2.486899575160351 × 10−13
250	$y_1$	2.462030579408747 × 10−12	7.105427357601002 × 10−15
	$y_2$	6.116351869422942 × 10−11	6.494360604847316 × 10−11
	$y_3$	1.620037437533028 × 10−12	0
500	$y_1$	2.613464999967619 × 10−13	6.661338147750939 × 10−16
	$y_2$	2.717115421546623 × 10−11	2.761169071163749 × 10−11
	$y_3$	1.421085471520200 × 10−13	1.705302565824240 × 10−13

As clearly seen in the numerical computations, even if taken step length or mesh point very small, due to the stiffness of the equations obtained from the formulation, one will still encounter problem of accuracy of results. This is exactly what happened in the results of both models II and IV presented in Tables 5 and 6, respectively. Especially, in the Table of values of model IV (Table 6), the two methods used are vying over accuracy. Though, the Ode45 from MATLAB has order higher than the method derived in (9) but both give the same accuracy of results. This account is for their RAS.

Table 5.

Absolute errors of numerical solutions of model II or Example 2.

$x$	$y_i$	Method (9) $|y(x)-y_n(x)|$	Method (11) $|y(x)-y_n(x)|$
5	$y_1$	1.838565004367581 × 10−00	1.838565000668211 × 10−01
	$y_2$	3.120693037317324 × 10−09	2.510347485440434 × 10−11
50	$y_1$	1.533692212025429 × 10−01	1.533692209989633 × 10−01
	$y_2$	1.703287466625625 × 10−08	2.785611741273897 × 10−10
150	$y_1$	2.343045476292816 × 10−01	2.343045475247055 × 10−01
	$y_2$	7.670507784496294 × 10−09	1.221849288413068 × 10−09
250	$y_1$	2.359885222398320 × 10−01	2.359885222103083 × 10−01
	$y_2$	1.665512172621675 × 10−11	2.526519438106334 × 10−09
350	$y_1$	2.249741755033404 × 10−01	2.249741754963486 × 10−01
	$y_2$	3.241382273699855 × 10−09	3.834998096863274 × 10−09
500	$y_1$	2.063157885989248 × 10−01	2.063157885982289 × 10−01
	$y_2$	5.508468348125462 × 10−09	5.566761274167220 × 10−09

Table 6.

Absolute errors of numerical solutions of model IV or Example 4.

$x$	$y_i$	Method (9) $|y(x)-y_n(x)|$	Ode 45 Method $|y(x)-y_n(x)|$
5	$y_1$	1.063274623902544 × 10−08	1.063274623902544 × 10−08
	$y_2$	4.380001476637586 × 10−09	4.380001476637586 × 10−09
	$y_3$	4.440000077855491 × 10−09	4.440000077855491 × 10−09
50	$y_1$	2.000011258829204 × 10−07	2.000011258829204 × 10−07
	$y_2$	4.380005454935290 × 10−09	4.380005454935263 × 10−09
	$y_3$	4.440000190123469 × 10−09	4.440000190123466 × 10−09
250	$y_1$	1.041638200582895 × 10−06	1.041638200582895 × 10−06
	$y_2$	4.380098253452721 × 10−09	4.380098253452721 × 10−09
	$y_3$	4.440009688640828 × 10−09	4.440009688640828 × 10−09
500	$y_1$	2.093684161152964 × 10−06	2.093684161152964 × 10−06
	$y_2$	4.380392587338717 × 10−09	4.380392587338717 × 10−09
	$y_3$	4.440038528554658 × 10−09	4.440038528554684 × 10−09

In order to obtain better results, the solutions of the problems of models III and IV of Figures 3 and 4 have been obtained. It is evident from the plots of the three Figures 8 to 10, that parts of the drugs, after absorption were transferred from the blood streams into tissues through the arterial blood stream from the tissues with fraction of drugs taken into the blood through the veins. These indicate the variation of drugs concentration in the arterial blood and venous blood. As shown from the plotted efficiency curves, it is clear that drug concentration in the arteries blood stream decreases rapidly and increases in the venous blood. This is a remarkable improvement to the compartmental modeling: in fact, once the model parameters have been evaluated for a certain administration, the model is able to predict the drug plasma concentration vary not only the dose, but also the infusion rate of the drug. In particular, the compartmental analysis has been extensively used and compared with the experimental data, taken after intravenous infusion ⁷ and bolus injection. ⁶

Figure 9.

Rapid process of drug delivery and concentration pattern within the body with different rate constants. (a) Graph of the Ode 45 Method. (b) Graph of Method (11).

Figure 10.

Rapid process of drug delivery and concentration pattern within the body with different rate constants. (a) Graph of Method(9). (b) Graph of Ode 45 Method (11).

Summary

In this section, we summarize our findings in the study especially the advantages of using block methods in the solution of stiff pharmacokinetic models. Such models are used to describe measurements, to categorize the pharmacological effect of different compounds, to simulate different dosing schedules (e.g. for first-in-human dose selection) and also to understand underlying mechanisms of disease and drug response. In the paper, we studied models of three types of problems on drug diffusion in the biological tissues. We considered the limitations arising from the diffusion processes taking place across the compartments in the general biological systems which can extensively be used for the drug diffusion found in pharmaceutical studies. Furthermore, the models generalized the properties of drug diffusion to flexibility of parameters and the importance of drug diffusion through different routes of drug administration. We are convinced, that mathematics has important impact on drug development and it is commonly believed that the role of mathematical modeling and simulations will further increase as efforts to establish individual clinical dosage plans and to standardize drug preparations have increased over the years due to some development from mathematical modeling. Important justification of these results are the comparison with curve fitting of the experimental data and the calculation of the areas under concentration curves or of the mean residence time of the drug in a given compartment or in the system as a whole.^6,7 Hence, we now conclude that from the numerical simulations performed on the stiff pharmacokinetic models in this paper on very large interval demonstrate the efficiency and exactness of the continuous implicit block hybrid one-step collocation methods compared to other numerical methods for blood distribution in the biological systems. The presented results will now widen the applicability of the continuous block implicit hybrid one-step collocation methods to diffusion process because of the good number of applications in the drug control, drug dosage and other related problems in pharmaceutical and biomedical industries.