A three-compartmental delay model is formulated to describe the pharmacokinetics of drugs subjected to both intravenous and oral doses with reabsorptions by the central compartment. Model dynamics are analyzed rigorously, and two equilibrium points are obtained to be locally asymptotically stable under certain conditions. Time delays used as lags in reabsorption of drugs by central compartment from other two compartments caused rebounds or peaks and fluctuations in the time profiles for amounts of drug in all the compartments. Sensitivity analysis revealed that elimination rates decrease the amounts in all compartments. Furthermore, reabsorption rates cause superimposition at the initial phases of the drug amount profiles; subsequently, the quantities decrease in compartment one and increase in compartments two and three, respectively.

1. Introduction

Pharmacokinetics (PK) is a branch of pharmacology (study of drug or medication action), derived from two Greek words pharmakon (drug) and kinetics (movement). In its simplest meaning, PK is the study of drug absorption, distribution, and elimination (related to excretion and metabolism) [1, 2]. Closely related to PK is pharmacodynamics (PD) which refers to the relationship between drug concentration at the site of action and the resulting effect, including the time course and intensity of therapeutic and adverse effects [1]. In the concept of many authors, there is no clear distinction between PK and PD. However, according to Holford and Sheiner [3], pharmacokinetics is described as “what the body does to the drug” and pharmacodynamics as “what the drug does to the body.” The most important contribution to drug research can be traced back to the pioneer work of Teorell [4] in his two famous articles which led to his title as “father of pharmacokinetics” [5]. In the articles, Toerell derived differential equations to describe the kinetics for distribution of drugs in the body and presented time-concentration curves graphically to depicts the relation.

Administrations of drug into body system are usually carried out through two most populous routes: enteral (oral, sublingual, and rectal) and parenteral (intravenous, intra-arterial injections, etc.) [6]. Entrohepatic circulation is one of the important drug transformation in the body and is a process in which administrated product excreted by liver into the bile is reabsorbed by the liver via portal blood stream [79]. The PK and PD of drugs after administration can be very complex as a result of various processes taking place that alter the amount of drug and/or concentrations in tissues and fluids [1]. Furthermore, measuring the amount of drug at the site of action and beyond is practically difficult if not impossible. Thus, there is a need for simplified body process representation in relation to the drug behaviour. One of the methods for body simplification process is the application of mathematical principles, using suitable compartmental models as introduced by Sheppard [10]. A compartment can represent a tissue, fluid, or group which is homogeneous and specific with respect to biological activities in the body with size determined by concentration or amount of material [11]. Based on the nature of PK of drugs in the body, compartments can be numbered and organs/tissues classified according to highly perfused and less well perfused. Highly perfused compartment is also called central compartment and comprised of the blood, heart, liver, lungs, and kidney or even the whole body. The less well perfused compartment, sometimes called peripheral compartment, includes fat tissue, muscle tissue, and skeletal cerebrospinal fluid [1]. Apparently, a compartmental model provides a comprehensive form of representing body tissues/organs fluids into one, two, or multicompartments, monitors the drug’s movement and behaviour, and predicts the time, amount, and concentration as it changes with time.

Generally, ordinary differential equations are among the most powerful mathematical tools used to describe models that are assumed instantaneous time transfer of material flow between and within compartments [1215]. However, for most biological systems, this inter/intraflow of materials does not occur immediately but with time delays [16, 17]. In particular, delay may cause dynamical instabilities that generate periodic solutions, and biologically, this may affect normal physiological functions [1821]. Many mathematical delay models of pharmacokinetics have been studied to give insight into the dynamics and pharmacodynamics of drugs. In PK models, delay can be used to represent time lag in transfer and absorption or reabsorption between one part of the body to another [79, 22]. For instance Steimer et al. in [8] considered two-compartmental model (body and gastrointestinal tract) with delay representing time lag for reabsorption of drug after biliary excretion from the body (first compartment) to gastrointestinal (second compartment). They found that their model exhibits qualitative agreement with the two basic experimental observations (rebounds and slow terminal kinetics) on PK of drugs undergoing enterohepatic circulation. Plusquellec and Bousquet [9] proposed another two-component model with time lag describing PK of drugs subject to enterohepatic circulation after an intravenous bolus, or an oral intake or a constant infusion. Using computer simulations, plasma level profiles that lead to rebounds and secondary peaks are obtained in accordance with experimental evidences. Labat et al. in [7] presented yet another two-component model with time delay for PK of drugs subject to enterohepatic circulation when reabsorption is repeated at unequal intervals. It was found that the first reabsorption peak area depends on time of administration and was noticeable only when a parameter was above a given threshold. In [22], analytic expression is presented for the solution of a delay differential system for any compartments with a single input and by convolution for all intakes. The result is applied to one-, two-, and three-compartmental models that provide insight into the PK of drug subjected to enterohepatic circulation among others.

From the literature presented above, it is evident that time delay play important roles in both the dynamics of PK models and physiological functions in the body. Motivated by these, in the current study, we extend three-compartmental models in [12, 23] by incorporating two delays, and , to represent time lags for reabsorption of drug at compartment 1 (central compartment) after elimination by compartments 2 and 3, respectively. Furthermore, our new PK delay model considers intravenous and oral doses of drug at the central and peripheral compartments, respectively, in addition to elimination in all the three compartments. Most of the parameter values used are obtained from published articles optimized from experimental data to validate the model. The objectives are to find expression that can determine the amount of drug at any given time with specified delay and investigate the effects of time lag, rate constants of transfer between/out of compartments. In what follows, the paper is organized as shown. In Section 2, the model formulation and description are given, followed by basic results in Section 3. Section 4 consists of equilibria and stability analysis, while numerical simulations and sensitivity analysis are presented in Section 5 to illustrate our equilibrium results and effects of variation of parameter values. Finally, the concluding remarks is presented in Section 6 followed by references.

2. Model Formulation and Description

The model consists of three compartments, with and representing the amounts of drug at time . We assume that compartment 1 (central compartment) represents highly perfused organs/tissues in the body such as the heart, liver, kidney, brain, or plasma. Compartment 2 less highly perfused organs like the gastrointestinal tract and bile while compartment 3 less perfused organs such as the adipose and skeletal muscle. The rate constants and are the kinetics constant of drug transfer between the central and compartment 2, while and are the kinetics constant of drug transfer between the central and compartment 3. There are drug eliminations (clearances) rates, , , and , respectively, from each compartment. We introduce time delays and , to represent the time lags in reabsorption of drugs in the central compartment from compartments 2 and 3, respectively. Thus drug leaving compartments 2 and 3 at time can be reabsorbed at central compartment at times and , respectively. Drug is introduced by intravenous injection at constant rate in the central compartment, and orally into compartments 2 and 3, respectively, at constant rates and . In all the cases of transport and elimination of drugs in the compartments, we assume first-order kinetics. The schematic diagram of the model in Figure 1, together with the model equations in the following systems describe the model dynamics.

The model (1)–(3) satisfy the nonnegative initial data defined as where belong to a Banach space of continuous functions, mapping the interval into () with .

Letting to be vector for the amount of drug in the body at time , the model systems (1), (2), and (3) can be written as where

3. Preliminary Results

Expressing the model equations in the form (5), it is well known from the fundamental theories of functional differential equations [24, 25] that model (1)–(3) admit a unique solution with the stated initial data.

Furthermore, since the model (1)–(3) represent biological process, we assume that all parameters and variables are positive at all times .

On the boundedness of solutions, we present the following result.

Theorem 1. Assuming the positivity of solutions with initial data, the solution set is bounded in the region

Proof. Assuming positivity of solution with initial data and letting be the total amount of drug in the body at time , adding the equations in models (1), (2), and (3), we have

4. Equilibria and Stability of the Model

Here, we present the equilibria solutions and their stabilities for the model. At equilibrium, we have , for . In the absence of intravenous and oral doses (), the trivial equilibrium exists and, from the model (1)–(3), is obtained to be

It is worth remarking here that the trivial equilibrium represent the initial amount of drug in the body parts when there is no intravenous and oral doses. With time, the existing amount of drug will be exhausted by the body metabolism and hence will be zero.

However, when there is introduction of either intravenous and/or oral doses (), the nontrivial equilibrium emerged and is obtained as follows. From (1), (2), and (3), at equilibrium, we get

Let . From (9),

Also from (13),

Substituting (15) and (17) in (11) and simplifying, gives

Therefore, from Equations (15)–(18), the unique equilibrium is given by

The nontrivial equilibrium always exists when . It is easy to see, based on the values of s, that this condition always holds true. We now determine an arbitrary transcendental equation for the model as follows, assuming a solution of the form , with an eigenvalue and a constant vector. From (5)

Taking as an identity vector, and for nontrivial solution, we have where is a identity matrix.

4.1. Stability of Trivial Equilibrium,

The stability of trivial equilibrium is stated here and proved.

Theorem 2. The trivial equilibrium is locally asymptotically stable for all values of delays and whenever , for all .

Proof. At , , and taking the determinant of Equation (21), we have When , Equation (22) becomes It can be seen that

This implies that all the eigenvalues of have negative real parts; hence, is stable for .

When , from (22), after simplification, we again have

Hence, has no positive real roots. Therefore, is absolutely stable for all delays and .

4.2. Stability of Nontrivial Equilibrium,

The equilibrium emerged whenever either or all of , for . With this condition, the determinant of transcendental Equation (21) is given as

Hence, Equation (26) can be express as where and are polynomials of degrees 3 and 2, respectively, defined as follows:

To discuss the stability of , we first study the distribution of roots for Equation (26). It can be seen from Equation (27) that

This indicates that Equation (27) has at least one positive root in . However, when both and , Equation (26) will be too complicated for analysis and any meaningful conclusion; thus, for further discussions, we will consider one case (, ), where is a stable interval for , using the approaches in [20, 26]. However, for the remaining two cases: case 2, , and case 3, , , where is a stable interval for , we will use numerical simulations to illustrate

Case 1: , .

In this case, we consider Equation (26) with in its stable interval and as the bifurcation parameter. Equation (26) will now be expressed as where and are defined as

Assuming that , is the root of (30). Substituting, simplifying, and separating the real and imaginary parts, we obtain

To get an expression for , we square both sides of (32) and (33) and add, so that

Without loss of generality, we assume that Equation (37) has three positive roots denoted as . For any given root, the critical delay is therefore given by

Therefore, when , for , then is a pair of purely imaginary roots of Equation (30) for . It is easy to see that the sequence is monotone increasing so that .

Let define

In order to be sure that are simple purely imaginary roots of Equation (30), we proceed by differentiating (21) with respect to to get

Substituting and simplifying, we have


Clearly, when is satisfied, it follows that .

According to the Hopf bifurcation theorem as stated in [26], this result can be summarized below.

Theorem 3. Suppose Equation (30) admits at least one positive root and is in its stable interval . The nontrivial equilibrium is locally asymptotically stable for any . Furthermore, the model (1)–(3) will undergo the Hopf bifurcation at when and family of periodic solutions bifurcate from the equilibrium when .

5. Numerical Simulations

In this section, we carry out some numerical experiments using parameter values as displayed in Table 1, obtained from experimental data in the literature to illustrate our theoretical results. In addition, we show the phenomenon of the Hopf bifurcation and sensitivity analysis of key parameters to illustrate effects of time delay for reabsorption of drugs, and eliminations in the model systems (1), (2), and (3).

Figure 2 illustrates the local asymptotic stability of trivial equilibrium using different delay values as proved in Theorem 2. It can be observed in both Figures 2(a) and 2(b) that when there is no intravenous and oral doses, irrespective of the delay values ( or and or ), the amount of drug in the compartments will be exhausted (tends to ) as time progresses.

5.1. Numerics for Stability of Equilibria

The dynamics of case 1, where and is illustrated in Figure 3 using time series and phase portrait graphs. In Figure 3(a), the local asymptotic stability of is shown when . However, when , loses its stability as can be seen in Figure 3(a), and family of periodic solutions bifurcates from the equilibrium point .

In Figure 4(a), the local stability of is depicted when case 2 is considered. When , , the nontrivial equilibrium , is shown to be locally asymptotically stable. The equilibrium then undergoes the Hopf bifurcation at ; thereafter, becomes unstable and family of solutions bifurcates from it as evident from Figure 4(b).

Similarly, case 3 is depicted in Figure 5. Here, the local asymptotic stability of the equilibrium is shown when in Figure 5(a), where are considered. However, when , loses its stability resulting in family of periodic solutions bifurcating from the equilibrium as displayed in Figure 5(b).

Remark 4. It is worth mentioning here that the practical implication of the stability for nontrivial equilibrium is that the amount of drugs will remain fixed at any given time in the three body parts considered whenever the reabsorption time from compartments 2 and 3 to 1 does not exceed certain critical value. However, when the absorption time exceeds such time, then the amounts of drug in the three body parts will be fluctuating.

5.2. Sensitivity Analysis of Parameters

In this section, we will show the variations of drug amounts as certain parameter values in the model (1)–(3) are varied using repeated computer simulations, so as to compare with other models in the literature.

We start with effect of time lag in reabsorption of drugs by compartment 1 (central compartment) from compartments 2 and 3. In Figures 6(a) and 6(b), we varied the values of with negligible in order to see the behaviour of drug’s amount as time progresses. It can be observed that as is increased, there are increase in number of rebounds or peaks in the time profiles of and in central and compartments 2, respectively. Similarly, when value is varied, as , with negligible , it can be seen from Figures 6(c) and 6(d) that the number of rebounds increased in the graphs of central and compartments 3 as time progresses. The implication of time lag on pharmacokinetics profile as reported in [8] is that the three phases can be distinguished. In the initial distribution phase, the amount of drug increases in the body with increase in time and delay, and the terminal distribution phase is where the amount stabilizes. Between these two phases, there is the third phase and transient phase, where the effect of delay is apparent with presence of rebounds or peaks. Therefore, as time delay is increased, as a results of the peaks, transient phase becomes longer. Finally, the whole pharmacokinetics process will be sensitive in its profile and duration as evident in Figures 6(a)6(d).

The effects of variation in reabsorption rates () and the drug’s amount profiles at the compartments are depicted in Figures 7(a)7(d). As can be seen, there are superimposition in all the initial phases of the profiles. However, as reabsorption rate from compartment 2 to 1, where is varied ( or ), the concentration of drug decreases in compartment 1 and increases in compartment 2 as shown in Figures 7(a) and 7(b), respectively. Similarly, when reabsorption rate from compartment 3 to 1, where is varied ( or ), the mass of drug decreases in compartment 1 while it increases in compartment 3 as evident from Figures 7(c) and 7(d), respectively.

The effects of elimination/excretion of drug from the three compartments are studied each by varying the corresponding rates. In Figure 8(a), the amount of drug profile in the central compartment is shown when the rate is varied () and the other two rates of eliminations are fixed as . In Figure 8(b) also, the amount profile in compartment 2 is displayed when the rate of elimination is varied () while . Similarly in Figure 8(c), the mass profile is shown when the elimination rate varies () or . In all the three cases, it can be observed that as the rates are increased, the respective amounts in the compartments decreases significantly as can be seen in Figures 8(a)8(c).

6. Concluding Remarks

In this article, a three-compartmental pharmacokinetics delay model is formulated and analyzed. The time delays represent the time lags in reabsorption of drug in the central compartment from two peripheral compartments. The model considers intravenous and oral administrations of drug to central and other two compartments, respectively. Furthermore, elimination of drug is considered in all the three compartments. Sensitivity analysis for effects of delays and parameter values obtained from published experimental data is performed numerically using MatLab dde23. Major findings of the research work are outlined below: (1)Two equilibrium solutions are obtained. The trivial equilibrium of the model, when there is no drug administered (intravenous or oral), is found to be locally asymptotically stable for all delay values. The biological meaning of this result is the initial amounts of drug in the compartments will vanish (approach zero) with time whenever there is no drug administered due to body metabolism. Peaks are also observed as reported in other models with delay. However, when there is either intravenous and/or oral doses, the nontrivial equilibrium emerged which is locally asymptotically stable for certain delay values. Thereafter, it becomes unstable and undergoes the Hopf bifurcation with family of periodic solutions bifurcating from the equilibrium. In this later case, the amount of drug can be determined within certain values of the delay as a result of delay effects in the model and thereafter, the amount fluctuates(2)Effects of time lags in reabsorption of drugs by the central compartment from compartments 2 and 3 as delays are varied are evident. It is observed that as delays are increased, there are increases in number of rebounds or peaks in the time profiles of amounts in all the compartments. These coincides with the pharmacokinetics profiles in published delay models. Thus, the pharmacokinetics process is sensitive in its profile and duration in the presence of time delay(3)When the reabsorption rates are increased, there are superimposition at the initial phases of the drug amount profiles; subsequently, the amount decreases in compartment 1 and increases in compartments 2 and 3, respectively(4)The effects of elimination/excretion of drug from the three compartments are studied each by varying the corresponding rates while other rates fixed. In all the compartments, it has been observed that as the rates are increased, the respective amount in the compartments decreases significantly.

Data Availability

Data are found in published articles, as stated in the manuscript.

Conflicts of Interest

The authors hereby declare that there is no conflicts of interest whatsoever.


A. S. Hassan acknowledges with thanks the support of the Research Support Directorate, University of South Africa, for the Research Grant Fellowship number 409000. He is also grateful for the authority of Bayero University Kano, Nigeria, for the one-year research leave.