# Effects of Drugs Reabsorption Time Delays for a Pharmacokinetics Model in Humans

Authors: Adamu Shitu Hassan; Justin M. W. Munganga
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

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.


---

## Body

## 1. Introduction

Pharmacokinetics (PK) is a branch of pharmacology (study of drug or medication action), derived from two Greek wordspharmakon (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 [7–9]. 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 [12–15]. 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 [18–21]. 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 [7–9, 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 n 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, τ1 and τ2, 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, withX1t,X2t, and X3t representing the amounts of drug at time t. 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 k12 and k21 are the kinetics constant of drug transfer between the central and compartment 2, while k13 and k31 are the kinetics constant of drug transfer between the central and compartment 3. There are drug eliminations (clearances) rates, k10, k20, and k30, respectively, from each compartment. We introduce time delays τ1 and τ2, 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 t can be reabsorbed at central compartment at times t+τ1 and t+τ2, respectively. Drug is introduced by intravenous injection at constant rate D1 in the central compartment, and orally into compartments 2 and 3, respectively, at constant rates D2 and D3. 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.
(1)dX1tdt=D1−k10+k12+k13X1t+k21X2t−τ1+k31X3t−τ2,(2)dX2tdt=D2−k21X2t−τ1−k20X2t+k12X1t,(3)dX3tdt=D3−k31X3t−τ2−k30X3t+k13X1t.Figure 1
Flow diagram of the pharmacokinetics delay model with reabsorption at the central compartment.The model (1)–(3) satisfy the nonnegative initial data defined as
(4)X1t=ψ1t,X2t=ψ2t,X3t=ψ3t,t∈−τ,0,where ψ1,ψ2,andψ3 belong to a Banach space of continuous functions, C mapping the interval −τ,0 into ℝ+3 (C−τ,0,ℝ+3) with τ=maxτ1,τ2.LettingXt=X1tX2tX3tT to be vector for the amount of drug in the body at time t, the model systems (1), (2), and (3) can be written as
(5)dXtdt=D+AXt+B1Xt−τ1+B2Xt−τ2, where
(6)D=D1000D2000D3,A=−K100k12−k200k130−k30,B1=0k2100−k210000,(7)B2=00k3100000−k31andXt=X1tX2tX3t,whereK1=k10+k12+k13.


## 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 X1t,X2t,X3t 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 t≥0.On the boundedness of solutions, we present the following result.Theorem 1.
Assuming the positivity of solutions with initial data, the solution setX1t,X2t,X3t is bounded in the region
(8)Ω=X1t,X2t,X3t∈ℝ+3∣XTt≤DKt,whereKt=mink10,k20,k30,D=D1+D2+D3.Proof.
Assuming positivity of solution with initial data and lettingXTt=X1t+X1t+X1t be the total amount of drug in the body at time t, adding the equations in models (1), (2), and (3), we have
(9)dXTtdt=D1+D2+D3−k10X1t−k20X2t−k30X3t,≤D−KtXTt,henceusingGronwallinequality,XTt≤DKtifXT0≤DKt.


## 4. Equilibria and Stability of the Model

Here, we present the equilibria solutions and their stabilities for the model. At equilibrium, we haveXit−τj=Xit=Xi∗, for i=1,2,3,j=1,2. In the absence of intravenous and oral doses (D1=D2=D3=0), the trivial equilibrium exists and, from the model (1)–(3), is obtained to be
(10)E0=X10,X20,X30=0,0,0.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 (D1≠D2≠D3≠0), the nontrivial equilibrium emerged and is obtained as follows. From (1), (2), and (3), at equilibrium, we get
(11)D1−K1X1∗+k21X2∗+k31X3∗=0,(12)D2−k21X2∗−k20X2∗+k12X1∗=0,(13)D3−k31X3∗−k30X3∗+k13X1∗=0.LetK2=k20+k21,K3=k30+k31. From (9),
(14)X2∗k20+k21=D2+k12X1∗,sothat(15)X2∗=D2+k12X1∗K2.Also from (13),
(16)X3∗k30+k31=D3+k13X1∗,hence(17)X3∗=D3+k13X1∗K3.Substituting (15) and (17) in (11) and simplifying, gives
(18)X1∗=D1+k21D2/K2+k31D3/K3K1−k21k12/K2−k31k13/K3.Therefore, from Equations (15)–(18), the unique equilibrium is given by
(19)E∗=X1∗,X2∗,X3∗.The nontrivial equilibrium always exists whenK1>k21k12/K2+k31k13/K3. It is easy to see, based on the values of kijs, that this condition always holds true. We now determine an arbitrary transcendental equation for the model as follows, assuming a solution of the form Xt=peλt, with λ an eigenvalue and p a constant vector. From (5)
(20)peλtλ=D+Apeλt+B1peλte−λτ1+B1peλte−λτ2.Takingp as an identity vector, and for nontrivial solution, we have
(21)λI−De−λt−A−B1e−λτ1−B1e−λτ2=0,where I is a 3×3 identity matrix.


### 4.1. Stability of Trivial Equilibrium,E0=0,0,0

The stability of trivial equilibrium is stated here and proved.Theorem 2.
The trivial equilibriumE0 is locally asymptotically stable for all values of delays τ1>0 and τ2>0 whenever Di=0, for all i=1,2,3.Proof.
AtE0=0,0,0, D1=D2=D3=0, and taking the determinant of Equation (21), we have
(22)gλ=λ+Kλ+k20+k21e−λτ1λ+k30+k31e−λτ2+k12−k21e−λτ1λ+k30+k31e−λτ2−k13k31e−λτ2λ+k20+k21e−λτ2=0,
Whenτ1=τ2=0, Equation (22) becomes
(23)Gλ=λ+Kλ+k20+k21λ+k30+k31+k12−k21λ+k30+k31−k13k31λ+k20+k21.
It can be seen that(24)G0=Kk20+k21k30+k31−k12k21k30+k31−k13k31k20+k21,=KK2K3−k12k21K3−k13k31K2=k10K2K3+k12k20K3+k13k30K2>0,andlimλ⟶∞Gλ=+∞.This implies that all the eigenvalues ofGλ=0 have negative real parts; hence, E0 is stable for τ1=τ2=0.Whenτ1>0,τ2>0, from (22), after simplification, we again have
(25)G0=k10K2K3+k12k20K3+k13k30K2>0,andlimλ⟶∞Gλ=+∞.Hence,Gλ=0 has no positive real roots. Therefore, E0 is absolutely stable for all delays τ1>0 and τ2>0.


### 4.2. Stability of Nontrivial Equilibrium,E∗=X1∗,X2∗,X3∗

The equilibriumE∗ emerged whenever either or all of Di>0, for i=1,2,3. With this condition, the determinant of transcendental Equation (21) is given as
(26)Hλ=λ−D1e−λt+K1λ−D2e−λt+k20+k21e−λτ1λ−D3e−λt+k30+k31e−λτ2−k12k21e−λτ1λ−D3e−λt+k30+k31e−λτ2−k13k31e−λτ2λ−D2e−λt+k20+k21e−λτ1=0.Hence, Equation (26) can be express as
(27)Hλ=Pλ+Qλ,where P and Q are polynomials of degrees 3 and 2, respectively, defined as follows:
(28)Pλ=λ3+p2λ2+p1λ+p0,Qλ=q2λ2+q1λ+q0,withp2=D3e−λt−k30+D2e−λt−k20+D1e−λt+K1,p1=−D2e−λt+k20−D3e−λt+k30−−D1e−λt+K−D3e−λt+k30−D2e−λt+k20,p0=−D1e−λt+K1−D3e−λt+k30−D2e−λt+k20,q2=−k21e−λτ1−k31e−λτ2,q1=k31e−λτ2−D2e−λt+k20+k21e−λτ1−D3e−λt+k30+k21k31e−λτ1e−λτ2−−D1e−λt+K1k21e−λτ1+k31e−λτ2−k12k21e−λτ1−k13k31e−λτ2,q0=−D1e−λt+K1D2e−λt−k20k31e−λτ2+−D1e−λt+K1−D3e−λt+k30k21e−λτ1+−D1e−λt+K1k21k31e−λτ1e−λτ2−k13k31e−λτ2D2e−λt−k20+k21e−λτ1.To discuss the stability ofE∗, we first study the distribution of roots for Equation (26). It can be seen from Equation (27) that
(29)H0=−D1D2D3+D1K2K3+K1K3D2+K1K2D3+k12k21k30+k12k21k31+k13k31k20+k13k31k21+D1D2K3+D1D3K2+k12k21D3+k13k31D2+K1D2D3+K1K2K3<0,andlimλ⟶∞Hλ=+∞.This indicates that Equation (27) has at least one positive root in 0,+∞. However, when both τ1>0 and τ2>0, Equation (26) will be too complicated for analysis and any meaningful conclusion; thus, for further discussions, we will consider one case (τ1>0, τ2∈0,τ20∗), where 0,τ20∗ is a stable interval for τ2, using the approaches in [20, 26]. However, for the remaining two cases: case 2, τ1=τ2=τ, and case 3, τ2>0, τ1∈0,τ10∗, where 0,τ20∗ is a stable interval for τ1, we will use numerical simulations to illustrateCase 1:τ1>0, τ2∈0,τ20∗.In this case, we consider Equation (26) with τ2 in its stable interval and τ1 as the bifurcation parameter. Equation (26) will now be expressed as
(30)Hλ=Pλ+Q1λe−λτ1+Q2λe−λτ2,where Q1 and Q2 are defined as
(31)Q1λ=−k21λ2+q3+q4λ+q5+q6,Q2λ=−k31λ2+q7λ+q8,withq3=k21−D3+k30−−D1+K1−k12,q4=q41e−λτ2,q41=k21k31,q5=−D1+K1−D3+k30k21,q6=q61e−λτ2,q61=−D1+K1k21k31−k13k31k21,q7=k31−D2+k20−−D1+K1k31−k13k31,q8=−D1+K1D2−k20k31+−k13k31−D2+k20.Assuming thatλ=iy, y∈ℝ+ is the root of (30). Substituting, simplifying, and separating the real and imaginary parts, we obtain
(32)P11cosy1τ1+Q11siny1τ1=N11,(33)Q11cosy1τ1+P11siny1τ1=M11,where(34)P11=y12k21+y1siny1τ2+q5+q61cosy1τ2,Q11=y12q3+y1q4cosy1τ2−q61siny1τ2,(35)M11=y13−y1p1+y12k31+q8siny1τ2−q7y1cosy1τ2,(36)N11=y12p2+p0+y12k31+q8cosy1τ2+q7y1siny1τ2.To get an expression fory1, we square both sides of (32) and (33) and add, so that
(37)y16+y14k212+p22+y12q3+q4cosy1τ22+∑i=14fiy1=0,where(38)f1=2y12k21y1siny1τ2+q5+q61cosy1τ2+q5+q61cosy1τ22,(39)f2=−2y1q3+q4cosy1τ2q61siny1τ2+q61siny1τ22,(40)f3=−2y13−y1p1+y12k31+q8siny1τ2−y1q7cosy1τ2+−y1p1+y12k31+q8siny1τ2−y1q7cosy1τ22,(41)f4=2y12p2y1k31y1+q7siny1τ2+p0+q8cosy1τ2+y1k31y1+q7siny1τ2+p0+q8cosy1τ22.Without loss of generality, we assume that Equation (37) has three positive roots denoted as y11,y12,andy13. For any given root, the critical delay is therefore given by
(42)τ1jm=1y1jcos−1Q11M11+P11N11P112+Q112+2πmy1j,j=1,2,3,m=0,1,2,.⋯Therefore, whenτ1=τ1jm, for j=1,2,3,m=0,1,2,⋯, then ±iy1j is a pair of purely imaginary roots of Equation (30) for τ2∈0,τ20∗. It is easy to see that the sequence τ1jm is monotone increasing so that limm⟶∞τ1jm=∞.Let define(43)τ10∗=min1≤j≤3,m≥0τ1jm.In order to be sure that±iy1j are simple purely imaginary roots of Equation (30), we proceed by differentiating (21) with respect to τ1 to get
(44)dλdτ1−1=3λ2+2p2λ+p1+e−λτ1−2k21λ+q3+q41e−λτ2−τ2e−λτ2q41+q61λe−λτ1−k21λ2+q3+q41e−λτ2λ+q5+q61e−λτ2+e−λτ2−2λk31+q7−τ2−λ2k31+λq7+q8λe−λτ1−k21λ2+q3+q41e−λτ2λ+q5+q61e−λτ2−τ1λ.Substitutingλ=iy1 and simplifying, we have
(45)dλdτ1−1=3λ2+2p2λ+p1+P11cosy1τ1+P12siny1τ1+2y1p2+P12cosy1τ1−P11siny1τ1iR11cosy1τ1+R12siny1τ1+R12cosy1τ1−R11siny1τ1i+Q11cosy1τ2+Q12siny1τ2+Q12cosy1τ2−Q11siny1τ2iR11cosy1τ1+R12siny1τ1+R12cosy1τ1−R11siny1τ1i−τ1iy1,whereP21=q3+q41cosy1τ2−τ2q41cosy1τ2−τ2q61cosy1τ2,P22=−2yk21−q41siny1τ2+τ2q41siny1τ2+τ2q61siny1τ2,Q21=q7−τ2y2k31−τ2q8,Q22=−2yk31−τ2yq7,R21=y2q3+y2q41cosy1τ2,R22=y3k21+y2q41siny1τ2.Hence(46)dReλdτ1τ1=τ10∗−1=Redλdτ1τ1=τ10∗−1=P3RP1R+P2R+Q3IQ1I+Q2IP3R2+Q3I2,whereP1R=−3y1∗+p1+P21cosy1∗τ1∗+P22siny1∗τ1∗,P2R=Q11cosy1∗τ2+Q22siny1∗τ2,P3R=R21cosy1∗τ1∗+R22siny1∗τ1∗,Q1I=2y1∗p2+P22cosy1∗τ1∗−P21siny1∗τ1∗,Q2I=Q22cosy1∗τ2−Q21siny1∗τ2.Clearly, whenP3RP1R+P2R+Q3IQ1I+Q2I≠0 is satisfied, it follows that Redλ/dτ1τ1=τ10∗−1≠0.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 τ2 is in its stable interval 0,τ20∗. The nontrivial equilibrium E∗ is locally asymptotically stable for any τ1∈0,τ10∗. Furthermore, the model (1)–(3) will undergo the Hopf bifurcation at E∗ when τ1=τ10∗ and family of periodic solutions bifurcate from the equilibrium when τ1>τ10∗.


## 4.1. Stability of Trivial Equilibrium,E0=0,0,0

The stability of trivial equilibrium is stated here and proved.Theorem 2.
The trivial equilibriumE0 is locally asymptotically stable for all values of delays τ1>0 and τ2>0 whenever Di=0, for all i=1,2,3.Proof.
AtE0=0,0,0, D1=D2=D3=0, and taking the determinant of Equation (21), we have
(22)gλ=λ+Kλ+k20+k21e−λτ1λ+k30+k31e−λτ2+k12−k21e−λτ1λ+k30+k31e−λτ2−k13k31e−λτ2λ+k20+k21e−λτ2=0,
Whenτ1=τ2=0, Equation (22) becomes
(23)Gλ=λ+Kλ+k20+k21λ+k30+k31+k12−k21λ+k30+k31−k13k31λ+k20+k21.
It can be seen that(24)G0=Kk20+k21k30+k31−k12k21k30+k31−k13k31k20+k21,=KK2K3−k12k21K3−k13k31K2=k10K2K3+k12k20K3+k13k30K2>0,andlimλ⟶∞Gλ=+∞.This implies that all the eigenvalues ofGλ=0 have negative real parts; hence, E0 is stable for τ1=τ2=0.Whenτ1>0,τ2>0, from (22), after simplification, we again have
(25)G0=k10K2K3+k12k20K3+k13k30K2>0,andlimλ⟶∞Gλ=+∞.Hence,Gλ=0 has no positive real roots. Therefore, E0 is absolutely stable for all delays τ1>0 and τ2>0.


## 4.2. Stability of Nontrivial Equilibrium,E∗=X1∗,X2∗,X3∗

The equilibriumE∗ emerged whenever either or all of Di>0, for i=1,2,3. With this condition, the determinant of transcendental Equation (21) is given as
(26)Hλ=λ−D1e−λt+K1λ−D2e−λt+k20+k21e−λτ1λ−D3e−λt+k30+k31e−λτ2−k12k21e−λτ1λ−D3e−λt+k30+k31e−λτ2−k13k31e−λτ2λ−D2e−λt+k20+k21e−λτ1=0.Hence, Equation (26) can be express as
(27)Hλ=Pλ+Qλ,where P and Q are polynomials of degrees 3 and 2, respectively, defined as follows:
(28)Pλ=λ3+p2λ2+p1λ+p0,Qλ=q2λ2+q1λ+q0,withp2=D3e−λt−k30+D2e−λt−k20+D1e−λt+K1,p1=−D2e−λt+k20−D3e−λt+k30−−D1e−λt+K−D3e−λt+k30−D2e−λt+k20,p0=−D1e−λt+K1−D3e−λt+k30−D2e−λt+k20,q2=−k21e−λτ1−k31e−λτ2,q1=k31e−λτ2−D2e−λt+k20+k21e−λτ1−D3e−λt+k30+k21k31e−λτ1e−λτ2−−D1e−λt+K1k21e−λτ1+k31e−λτ2−k12k21e−λτ1−k13k31e−λτ2,q0=−D1e−λt+K1D2e−λt−k20k31e−λτ2+−D1e−λt+K1−D3e−λt+k30k21e−λτ1+−D1e−λt+K1k21k31e−λτ1e−λτ2−k13k31e−λτ2D2e−λt−k20+k21e−λτ1.To discuss the stability ofE∗, we first study the distribution of roots for Equation (26). It can be seen from Equation (27) that
(29)H0=−D1D2D3+D1K2K3+K1K3D2+K1K2D3+k12k21k30+k12k21k31+k13k31k20+k13k31k21+D1D2K3+D1D3K2+k12k21D3+k13k31D2+K1D2D3+K1K2K3<0,andlimλ⟶∞Hλ=+∞.This indicates that Equation (27) has at least one positive root in 0,+∞. However, when both τ1>0 and τ2>0, Equation (26) will be too complicated for analysis and any meaningful conclusion; thus, for further discussions, we will consider one case (τ1>0, τ2∈0,τ20∗), where 0,τ20∗ is a stable interval for τ2, using the approaches in [20, 26]. However, for the remaining two cases: case 2, τ1=τ2=τ, and case 3, τ2>0, τ1∈0,τ10∗, where 0,τ20∗ is a stable interval for τ1, we will use numerical simulations to illustrateCase 1:τ1>0, τ2∈0,τ20∗.In this case, we consider Equation (26) with τ2 in its stable interval and τ1 as the bifurcation parameter. Equation (26) will now be expressed as
(30)Hλ=Pλ+Q1λe−λτ1+Q2λe−λτ2,where Q1 and Q2 are defined as
(31)Q1λ=−k21λ2+q3+q4λ+q5+q6,Q2λ=−k31λ2+q7λ+q8,withq3=k21−D3+k30−−D1+K1−k12,q4=q41e−λτ2,q41=k21k31,q5=−D1+K1−D3+k30k21,q6=q61e−λτ2,q61=−D1+K1k21k31−k13k31k21,q7=k31−D2+k20−−D1+K1k31−k13k31,q8=−D1+K1D2−k20k31+−k13k31−D2+k20.Assuming thatλ=iy, y∈ℝ+ is the root of (30). Substituting, simplifying, and separating the real and imaginary parts, we obtain
(32)P11cosy1τ1+Q11siny1τ1=N11,(33)Q11cosy1τ1+P11siny1τ1=M11,where(34)P11=y12k21+y1siny1τ2+q5+q61cosy1τ2,Q11=y12q3+y1q4cosy1τ2−q61siny1τ2,(35)M11=y13−y1p1+y12k31+q8siny1τ2−q7y1cosy1τ2,(36)N11=y12p2+p0+y12k31+q8cosy1τ2+q7y1siny1τ2.To get an expression fory1, we square both sides of (32) and (33) and add, so that
(37)y16+y14k212+p22+y12q3+q4cosy1τ22+∑i=14fiy1=0,where(38)f1=2y12k21y1siny1τ2+q5+q61cosy1τ2+q5+q61cosy1τ22,(39)f2=−2y1q3+q4cosy1τ2q61siny1τ2+q61siny1τ22,(40)f3=−2y13−y1p1+y12k31+q8siny1τ2−y1q7cosy1τ2+−y1p1+y12k31+q8siny1τ2−y1q7cosy1τ22,(41)f4=2y12p2y1k31y1+q7siny1τ2+p0+q8cosy1τ2+y1k31y1+q7siny1τ2+p0+q8cosy1τ22.Without loss of generality, we assume that Equation (37) has three positive roots denoted as y11,y12,andy13. For any given root, the critical delay is therefore given by
(42)τ1jm=1y1jcos−1Q11M11+P11N11P112+Q112+2πmy1j,j=1,2,3,m=0,1,2,.⋯Therefore, whenτ1=τ1jm, for j=1,2,3,m=0,1,2,⋯, then ±iy1j is a pair of purely imaginary roots of Equation (30) for τ2∈0,τ20∗. It is easy to see that the sequence τ1jm is monotone increasing so that limm⟶∞τ1jm=∞.Let define(43)τ10∗=min1≤j≤3,m≥0τ1jm.In order to be sure that±iy1j are simple purely imaginary roots of Equation (30), we proceed by differentiating (21) with respect to τ1 to get
(44)dλdτ1−1=3λ2+2p2λ+p1+e−λτ1−2k21λ+q3+q41e−λτ2−τ2e−λτ2q41+q61λe−λτ1−k21λ2+q3+q41e−λτ2λ+q5+q61e−λτ2+e−λτ2−2λk31+q7−τ2−λ2k31+λq7+q8λe−λτ1−k21λ2+q3+q41e−λτ2λ+q5+q61e−λτ2−τ1λ.Substitutingλ=iy1 and simplifying, we have
(45)dλdτ1−1=3λ2+2p2λ+p1+P11cosy1τ1+P12siny1τ1+2y1p2+P12cosy1τ1−P11siny1τ1iR11cosy1τ1+R12siny1τ1+R12cosy1τ1−R11siny1τ1i+Q11cosy1τ2+Q12siny1τ2+Q12cosy1τ2−Q11siny1τ2iR11cosy1τ1+R12siny1τ1+R12cosy1τ1−R11siny1τ1i−τ1iy1,whereP21=q3+q41cosy1τ2−τ2q41cosy1τ2−τ2q61cosy1τ2,P22=−2yk21−q41siny1τ2+τ2q41siny1τ2+τ2q61siny1τ2,Q21=q7−τ2y2k31−τ2q8,Q22=−2yk31−τ2yq7,R21=y2q3+y2q41cosy1τ2,R22=y3k21+y2q41siny1τ2.Hence(46)dReλdτ1τ1=τ10∗−1=Redλdτ1τ1=τ10∗−1=P3RP1R+P2R+Q3IQ1I+Q2IP3R2+Q3I2,whereP1R=−3y1∗+p1+P21cosy1∗τ1∗+P22siny1∗τ1∗,P2R=Q11cosy1∗τ2+Q22siny1∗τ2,P3R=R21cosy1∗τ1∗+R22siny1∗τ1∗,Q1I=2y1∗p2+P22cosy1∗τ1∗−P21siny1∗τ1∗,Q2I=Q22cosy1∗τ2−Q21siny1∗τ2.Clearly, whenP3RP1R+P2R+Q3IQ1I+Q2I≠0 is satisfied, it follows that Redλ/dτ1τ1=τ10∗−1≠0.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 τ2 is in its stable interval 0,τ20∗. The nontrivial equilibrium E∗ is locally asymptotically stable for any τ1∈0,τ10∗. Furthermore, the model (1)–(3) will undergo the Hopf bifurcation at E∗ when τ1=τ10∗ and family of periodic solutions bifurcate from the equilibrium when τ1>τ10∗.


## 5. Numerical Simulations

In this section, we carry out some numerical experiments using parameter values as displayed in Table1, 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).Table 1
Baseline values for the parameters of model (1)–(3).
ParameterDescriptionValueReferenceD1Rate of intravenous dose toX150mgmin−1[25]D2Rate of oral dose toX240mgmin−1[26]D3Rate of oral dose toX340mgmin−1Assumedk10Rate of central body elimination of drug0.172min−1[2]k12Rate of transporting drug from central to compartment 20.373min−1[2]k13Rate of transporting drug from central to compartment 30.0367min−1[2]k21Reabsorption rate of drug from compartment 1 to central0.103min−1[2]k20Rate of drug elimination from compartment 20.00828min−1[2]k30Rate of drug elimination from compartment 30.00784min−1[2]k31Reabsorption rate of drug from compartment 3 to central0.0124min−1[2]τ1Time lag for reabsorption of drug from compartment 2 to 1[1,50] min[11]τ2Time lag for reabsorption of drug from compartment 3 to 1[1,90] min[11]Figure2 illustrates the local asymptotic stability of trivial equilibrium E0 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 (τ1=5 or 100 and τ2=4 or 100), the amount of drug in the compartments will be exhausted (tends to 0) as time progresses.Figure 2
Numerical simulations displaying the local stability of trivial equilibriumE0, using parameter values in Table 1, using (a) small delays (τ1=5,τ2=4) and (b) large delays (τ1=100,τ2=100).
(a)(b)


### 5.1. Numerics for Stability of Equilibria

The dynamics of case 1, whereτ1>0 and τ2∈0,τ20∗ is illustrated in Figure 3 using time series and phase portrait graphs. In Figure 3(a), the local asymptotic stability of E∗ is shown when τ1=35<τ10∗=44.50,τ2=20. However, when τ1=45.60>τ10∗=44.50,τ2=20, E∗ loses its stability as can be seen in Figure 3(a), and family of periodic solutions bifurcates from the equilibrium point E∗.Figure 3
Numerical simulations for the dynamics of nontrivial equilibriumE∗ showing case 1 with time series and phase portrait graphs, in (a) local stability when τ1=35<τ10∗=44.50,τ2=15 and (b) instability for τ1=45.60>τ10∗=44.50,τ2=20, showing the Hopf bifurcation, using parameter values in Table 1.
(a)(b)In Figure4(a), the local stability of E∗ is depicted when case 2 is considered. When τ2=29<τ20∗=33.23, τ1=40∈0,τ10∗, the nontrivial equilibrium E∗, is shown to be locally asymptotically stable. The equilibrium then undergoes the Hopf bifurcation at τ2=τ20∗=33.23; thereafter, E∗ becomes unstable and family of solutions bifurcates from it as evident from Figure 4(b).Figure 4
Numerical simulations for the local stability of trivial equilibriumE∗, displaying case 2, using parameter values in Table 1: (a) τ1=40,τ2=29, (b) τ1=43,τ2=34.80.
(a)(b)Similarly, case 3 is depicted in Figure5. Here, the local asymptotic stability of the equilibrium E∗ is shown when in Figure 5(a), where τ1=τ2=25<τc∗=42.35 are considered. However, when τ=44.35>τc∗, E∗ 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 equilibriumE∗ 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.Figure 5
Numerical simulations for the local stability of trivial equilibriumE∗, displaying case 3, using parameter values in Table 1: (a) τ1=τ2=25, (b) τ1=τ2=44.35.
(a)(b)


### 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 Figures6(a) and 6(b), we varied the values of τ1=12,27,40 with negligible τ2=0.1 in order to see the behaviour of drug’s amount as time progresses. It can be observed that as τ1 is increased, there are increase in number of rebounds or peaks in the time profiles of X1 and X2 in central and compartments 2, respectively. Similarly, when value τ2 is varied, as τ2=30,60,90, with negligible τ1=0.1, 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).Figure 6
Time series for effects of time delay in reabsorption of drug from compartments 2 and 3 to the central compartment using parameter values in Table1, except for (a) and (b) τ1=12,27,40 and τ2=0.1 (c) and (d) τ2=30,60,90 and τ1=0.1.
(a)(b)(c)(d)The effects of variation in reabsorption rates (k21,k31) 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 k21 is varied (k21=0.0103,0.0206,0.1030 or k31=0), 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 k31 is varied (k31=0.00124,0.00248,0.0124 or k21=0), the mass of drug decreases in compartment 1 while it increases in compartment 3 as evident from Figures 7(c) and 7(d), respectively.Figure 7
Time series for effect of reabsorption rates using parameter values in Table1, except, τ1=τ2=30, in (a) and (b) k21=0.0103,0.0206,0.1030 and k31=0. In (c) and (d) k31=0.00124,0.00248,0.0124 and k21=0.
(a)(b)(c)(d)The effects of elimination/excretion of drug from the three compartments are studied each by varying the corresponding rates. In Figure8(a), the amount of drug profile in the central compartment is shown when the rate k10 is varied (k10=0.0208,0.0416,0.208) and the other two rates of eliminations are fixed as k20=k30=0. In Figure 8(b) also, the amount profile in compartment 2 is displayed when the rate of elimination k20 is varied (k20=0.000828,0.001656,0.00828) while k10=k30=0. Similarly in Figure 8(c), the mass profile is shown when the elimination rate k30 varies (k30=0.000784,0.001568,0.00784) or k10=k20=0. 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).Figure 8
Time series for effects of elimination rates of drug using parameter values in Table1, except τ1=τ2=30, in (a) k10=0.0208,0.0416,0.208 and k20=k30=0, (b) k20=0.000828,0.001656,0.00828 and k10=k30=0, (c) k30=0.000784,0.001568,0.00784 and k10=k20=0.
(a)(b)(c)


## 5.1. Numerics for Stability of Equilibria

The dynamics of case 1, whereτ1>0 and τ2∈0,τ20∗ is illustrated in Figure 3 using time series and phase portrait graphs. In Figure 3(a), the local asymptotic stability of E∗ is shown when τ1=35<τ10∗=44.50,τ2=20. However, when τ1=45.60>τ10∗=44.50,τ2=20, E∗ loses its stability as can be seen in Figure 3(a), and family of periodic solutions bifurcates from the equilibrium point E∗.Figure 3
Numerical simulations for the dynamics of nontrivial equilibriumE∗ showing case 1 with time series and phase portrait graphs, in (a) local stability when τ1=35<τ10∗=44.50,τ2=15 and (b) instability for τ1=45.60>τ10∗=44.50,τ2=20, showing the Hopf bifurcation, using parameter values in Table 1.
(a)(b)In Figure4(a), the local stability of E∗ is depicted when case 2 is considered. When τ2=29<τ20∗=33.23, τ1=40∈0,τ10∗, the nontrivial equilibrium E∗, is shown to be locally asymptotically stable. The equilibrium then undergoes the Hopf bifurcation at τ2=τ20∗=33.23; thereafter, E∗ becomes unstable and family of solutions bifurcates from it as evident from Figure 4(b).Figure 4
Numerical simulations for the local stability of trivial equilibriumE∗, displaying case 2, using parameter values in Table 1: (a) τ1=40,τ2=29, (b) τ1=43,τ2=34.80.
(a)(b)Similarly, case 3 is depicted in Figure5. Here, the local asymptotic stability of the equilibrium E∗ is shown when in Figure 5(a), where τ1=τ2=25<τc∗=42.35 are considered. However, when τ=44.35>τc∗, E∗ 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 equilibriumE∗ 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.Figure 5
Numerical simulations for the local stability of trivial equilibriumE∗, displaying case 3, using parameter values in Table 1: (a) τ1=τ2=25, (b) τ1=τ2=44.35.
(a)(b)


## 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 Figures6(a) and 6(b), we varied the values of τ1=12,27,40 with negligible τ2=0.1 in order to see the behaviour of drug’s amount as time progresses. It can be observed that as τ1 is increased, there are increase in number of rebounds or peaks in the time profiles of X1 and X2 in central and compartments 2, respectively. Similarly, when value τ2 is varied, as τ2=30,60,90, with negligible τ1=0.1, 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).Figure 6
Time series for effects of time delay in reabsorption of drug from compartments 2 and 3 to the central compartment using parameter values in Table1, except for (a) and (b) τ1=12,27,40 and τ2=0.1 (c) and (d) τ2=30,60,90 and τ1=0.1.
(a)(b)(c)(d)The effects of variation in reabsorption rates (k21,k31) 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 k21 is varied (k21=0.0103,0.0206,0.1030 or k31=0), 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 k31 is varied (k31=0.00124,0.00248,0.0124 or k21=0), the mass of drug decreases in compartment 1 while it increases in compartment 3 as evident from Figures 7(c) and 7(d), respectively.Figure 7
Time series for effect of reabsorption rates using parameter values in Table1, except, τ1=τ2=30, in (a) and (b) k21=0.0103,0.0206,0.1030 and k31=0. In (c) and (d) k31=0.00124,0.00248,0.0124 and k21=0.
(a)(b)(c)(d)The effects of elimination/excretion of drug from the three compartments are studied each by varying the corresponding rates. In Figure8(a), the amount of drug profile in the central compartment is shown when the rate k10 is varied (k10=0.0208,0.0416,0.208) and the other two rates of eliminations are fixed as k20=k30=0. In Figure 8(b) also, the amount profile in compartment 2 is displayed when the rate of elimination k20 is varied (k20=0.000828,0.001656,0.00828) while k10=k30=0. Similarly in Figure 8(c), the mass profile is shown when the elimination rate k30 varies (k30=0.000784,0.001568,0.00784) or k10=k20=0. 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).Figure 8
Time series for effects of elimination rates of drug using parameter values in Table1, except τ1=τ2=30, in (a) k10=0.0208,0.0416,0.208 and k20=k30=0, (b) k20=0.000828,0.001656,0.00828 and k10=k30=0, (c) k30=0.000784,0.001568,0.00784 and k10=k20=0.
(a)(b)(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 MatLabdde23. 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.



---

*Source: 3439887-2021-12-26.xml*