Journal List > Transl Clin Pharmacol > v.27(2) > 1138101

D'Argenio and Bae: Analytical solution of linear multi-compartment models with non-zero initial condition and its implementation with R

Abstract

The analytical solution for multi-compartment models with a non-zero initial condition is complex because of the inter-compartmental transfer. An elegant solution and its implementation in the ‘wnl' R package can be useful in solving examples of textbooks and developing software of therapeutic drug monitoring, pharmacokinetic simulation, and parameter estimation. This solution uses Laplace transformation, convolution, matrix inversion, and the fact that the general solution of an inhomogeneous ordinary differential equation is the sum of a homogenous and a particular solution, together.

An analytical solution is not a numerical solution, and it can calculate concentrations at any time point without calculating the intervening time points (if the input function was not changed meanwhile). On the contrary, it is necessary to calculate concentrations of intervening many time points if one has only a numerical solution (in other words, if you use numerical integration).
If a drug is administered without previous dosing and not an endogenous substance, the drug concentration at time 0 (initial condition) will be zero. However, if the drug is cleared insufficiently from the previous dose or an endogenous substance, the drug concentrations of each compartment at time 0 (the initial condition) would not be zero. This condition is called non-zero initial condition here.
With zero initial values, it is not so difficult to have an analytical solution. And, if the pharmacokinetic model is a one-compartment model, Non-zero initial value problem can be solved without much difficulty using the rule of superposition for the drug amount (or concentration).
However, if the pharmacokinetic model is a multi-compartment model and the initial condition is not-zero for all compartments, one can not simply use the rule of superposition for the concentrations because the drug of each compartment moves between compartments.
For the simulation, numerical integration is often sufficient. However, for the fitting or estimation, an analytical solution is usually better because of the speed. Solutions with a single dose or zero-initial value are easily found in textbooks. However, solutions for the multi-compartment model with non-zero initial value are harder to find. Authors would like to present the most elegant way to solve this kind of problems analytically.

Case presentation

[Dosing history]
  • At 0-hour, 100 mg of a drug was administered by intravenous bolus.

  • At 24-hour, 150 mg of a drug was infused intravenously with the rate of 50 mg/hr.

  • At 48-hour, 100 mg of a drug was administered orally.

[Measurement time points]
  • The concentrations of a drug at 0, 1, 2, 4, 8, 12 hours after each dose with given pharmacokinetic parameters.

[Exercise 1. One-compartment model]
  • If the drug is disposed by a one-compartment model with the pharmacokinetic parameters of Ka = 1, Ke = 0.1, F = 1, and V = 1, what will be the concentrations at the measurement time points?

[Exercise 2. Two-compartment model]
  • If the drug is disposed by a two-compartment model with the pharmacokinetic parameters of Ka = 1, Ke = K10 = 0.1, K12 = 3, K21 = 1, F = 1, and V = 1, what will be the concentrations at the measurement time points?

[Exercise 3. Three-compartment model]
  • If the drug is disposed by a three-compartment model with the pharmacokinetic parameters of Ka = 1, Ke = K10 = 0.1, K12 = 3, K21 = 1, K13 = 2, K31 = 0.5, F = 1, and V = 1, what will be the concentrations at the measurement time points?

The easiest way to solve the above problem is using the ‘wnl' R package. This package can be used to solve all of the pharmacokinetic and pharmacodynamic exercises in the pharmacokinetic textbooks or manuals. This can also be called in other packages for therapeutic drug monitoring, parameter estimation, pharmacokinetic or pharmacodynamic simulation. This package is for deterministic simulations, but not for stochastic ones. For a stochastic simulation, one needs to generate random pharmacokinetic parameters using multivariate distributions and random residuals.
To install this from the default repository, use the following command.
install.packages(“wnl”)
require(wnl)
  • # Data preparation

  • DAT

  • DAT2 = ExpandDH(DAT) ; DAT2

  • # Solution 1

  • X1 = Comp1(Ke = 0.1, Ka = 1, DAT2) ; X1

  • matplot(DAT2[, “TIME”], X1, type = “l”)

  • # Solution 2

  • Sol = SolComp2(K10 = 0.1, K12 = 3, K21 = 1)

  • X2 = nComp(Sol, Ka = 1, DAT2) ; X2

  • matplot(DAT2[, “TIME”], X2, type = “l”)

  • # Solution 3

  • Sol = SolComp3(K10 = 0.1, K12 = 3, K21 = 1, K13 = 2, K31 = 0.5)

  • X3 = nComp(Sol, Ka = 1, DAT2) ; X3

  • matplot(DAT2[, “TIME”], X3, type = “l”)

The results of the above R script were confirmed with ADAPT 5 and NONMEM 7.4.3.
One can see all the source script with the following commands, and all the functions have only 20 to 40 lines for each.
  • ExpandDH # Expand Dosing History to include all non-differentiable points

  • Comp1 # One-compartment model with pharmacokinetic parameters

  • nComp # Multi-compartment model with lambdas and coefficients

  • SolComp2 # Solver for lambdas and coefficients of a 2-compartment model

  • SolComp3 # Solver for lambdas and coefficients of a 3-compartment model

There are many ways to solve this kind of problems. The strategy adopted here is that each amount at time ti will be calculated using the elapsed time (delta t) from the previous time point ti-1 and an initial condition equal to the condition (amounts of each compartment) of the previous time ti-1 point. This strategy is used for software such as NONMEM and ADAPT.
The list of terms related to the differential equation but not explained here are ordinary differential equation, initial value problem, coupled linear differential equation, homogenous differential equation, inhomogeneous differential equation, homogenous solution, particular solution, impulse response function, Laplace transformation, and convolution. The definitions and meaning of these terms can be found in textbooks or on the internet without difficulty.

One-compartment model with non-zero initial condition

After intravenous bolus injection, the drug amount at time t, X(t) in the central compartment (or concentration if you divide the amount by the volume of distribution) can be calculated using the following equation.
eq. 1
Xt=X(0)eKettcp-27-43-e001.jpg
where X(0) is the initial condition, the amount in the body at the time 0, and Ke the elimination rate constant.
Here e−Ket is called the impulse response function, or the response function to the unit impulse.
Eq. 1 is the solution to the following differential equation (Eq. 2).
Eq. 2
dXdt=KeXtcp-27-43-e002.jpg
If the drug is infused intravenously with the rate R (Eq. 3), the solution is Eq. 4.
Eq. 3
dXdt=KeX+Rtcp-27-43-e003.jpg
Eq. 4
Xt=X0eKet+RKe1eKettcp-27-43-e004.jpg
Eq. 4 is also viewed as a homogenous solution + a particular solution because Eq. 2 is a homogenous differential equation and Eq. 3 is an inhomogeneous differential equation.
The followings are true in general.
Eq. 5
dXdt=AX+f(t)withinitialconditionXotcp-27-43-e005.jpg
Eq. 6
X(t)=XH(t)+XP(t)tcp-27-43-e006.jpg
Eq. 7
X(t)=expAtXo+expAt*f(t)tcp-27-43-e007.jpg
Eq. 8
X(t)=expAtXo+0texpAtτf(τ)dτtcp-27-43-e008.jpg
where
A: coefficient matrix
X(t)tcp-27-43-e042.jpg: general solution of the differential equation dXdt=AX+f(t)tcp-27-43-e043.jpg
XH(t)=expAtXotcp-27-43-e044.jpg: solution of the homogenous differential equation dXdt=AXtcp-27-43-e045.jpg
exp(At): impulse response function
f(t)tcp-27-43-e046.jpg: input function
XP(t)tcp-27-43-e047.jpg: particular solution of the inhomogeneous differential equation
*: convolution operator
If a drug is infused and administered orally, the differential equation is
Eq. 9
dXdt=KeX+R+KaXg(0)eKattcp-27-43-e009.jpg
where Ka is the absorption rate constant, and Xg(0) is the orally administered dose multiplied by its bioavailability.
Then the solution is
Eq. 10
X(t)=XHt+XP,IVt+XP,POt=expKetX(0)+0texp[Ketτ]Rdτ+0texp[Ketτ]KaXg(0)eKaτdτtcp-27-43-e010.jpg
Eq. 11
Xt=X(0)eKet+RKe1eKet+KaXg(0)KaKeeKeteKattcp-27-43-e011.jpg
This is interestingly the summation of each solution for intravenous bolus, infusion, and oral dosing.
If Ka = Ke = K, the solution becomes
Eq. 12
Xt=X(0)eKt+RK1eKt+Xg(0)KteKttcp-27-43-e012.jpg

Two-compartment model with non-zero initial condition

The differential equation for a typical two-compartment model is
Eq. 13
dXdt=(Ke+Kcp)KpcKcpKpcX+R0+KaXg(0)eKat0tcp-27-43-e013.jpg
where
Ke = K10: the elimination rate constant from the central compartment
Kcp = K12: the rate constant from the central to the peripheral compartment
Kpc = K21: the rate constant from the peripheral to the central compartment
R: infusion rate to the central compartment. Infusion is given to the central compartment only.
Ka: the absorption rate constant from the gut compartment
Xg(0): amount in the gut compartment at time 0. It is after the multiplication with bioavailability.
The solution of the above would be the following.
Eq. 14
X=XHt+XP,IVt+XP,POt=expAtX(0)+0texp[Atτ]R0dτ+0texp[Atτ]KaXg(0)eKaτ0dτwhenA=(Ke+Kcp)KpcKcpKpctcp-27-43-e014.jpg
Eq. 15
expAt=L1(SIA)1tcp-27-43-e015.jpg
L−1 indicates inverse Laplace transformation.
For the matrix within the brackets,
Eq. 16
SIA=S+(Ke+Kcp)KpcKcpS+Kpctcp-27-43-e016.jpg
is inverted to have
Eq. 17
SIA1=S+KpcKpcKcpS+Ke+KcpS2+Ke+Kcp+KpcS+KeKpctcp-27-43-e017.jpg
The denominator is a quadratic equation, and if the roots are −λ1 and −λ2,
Eq. 18
SIA1=S+KpcS+λ1S+λ2KpcS+λ1S+λ2KcpS+λ1S+λ2S+Ke+KcpS+λ1S+λ2tcp-27-43-e018.jpg
where
Eq. 19
λ1=(Ke+ Kcp+Kpc)+(Ke+Kcp+Kpc)24KeKpc2tcp-27-43-e019.jpg
Eq. 20
λ2=Ke+ Kcp+Kpc(Ke+Kcp+Kpc)24KeKpc2tcp-27-43-e020.jpg
Here, we present only the case that roots are two distinct real values, because this is practically true in real situation.
Eq. 21
SIA1=C1S+λ1+C2S+λ2C3S+λ1+C4S+λ2C5S+λ1+C6S+λ2C7S+λ1+C8S+λ2tcp-27-43-e021.jpg
where
C1=λ1Kpcλ1λ2tcp-27-43-e048.jpg C2=Kpcλ2λ1λ2tcp-27-43-e049.jpg C3=Kpcλ1λ2tcp-27-43-e050.jpg C4=Kpcλ1λ2tcp-27-43-e051.jpg
C5=Kcpλ1λ2tcp-27-43-e052.jpg C6=Kcpλ1λ2tcp-27-43-e053.jpg C7=λ1KeKcpλ1λ2tcp-27-43-e054.jpg C8=Ke+Kcpλ2λ1λ2tcp-27-43-e055.jpg
Therefore,
Eq. 22
expAt=λ1Kpcλ1λ2eλ1t+Kpcλ2λ1λ2eλ2tKpcλ1λ2eλ1t+Kpcλ1λ2eλ2tKcpλ1λ2eλ1t+Kcpλ1λ2eλ2tλ1KeKcpλ1λ2eλ1t+Ke+Kcpλ2λ1λ2eλ2ttcp-27-43-e022.jpg
Eq.23
expAt=1λ1λ2(λ1Kpc)eλ1t+(Kpcλ2)eλ2tKpc·eλ1t+Kpc·eλ2tKcp·eλ1t+Kcp·eλ2t(λ1KeKcp)eλ1t+(Ke+Kcpλ2)eλ2ttcp-27-43-e023.jpg
Therefore, the homogenous solution is
Eq. 24
XHt=expAtX(0)tcp-27-43-e024.jpg
Then the particular solution for intravenous dosing is
Eq. 25
XP,IVt=Rλ1λ20t(λ1Kpc)eλ1(tτ)+(Kpcλ2)eλ2(tτ)Kcp·eλ1(tτ)+Kcp·eλ2(tτ)dτtcp-27-43-e025.jpg
Because the infusion rate for the peripheral compartment is zero, we can only consider the first column of exp(At).
If we integrate Eq. 25,
Eq. 26
XP,IVt=Rλ1λ2λ1Kpcλ1eλ1(tτ)+Kpcλ2λ2eλ2(tτ)Kcpλ1eλ1(tτ)+Kcpλ2eλ2(tτ)0ttcp-27-43-e026.jpg
and
Eq. 27
XP, IVt=Rλ1λ2λ1Kpcλ1(1eλ1t)+Kpcλ2λ2(1eλ2t)Kcpλ1(1eλ1t)+Kcpλ2(1eλ2t)tcp-27-43-e027.jpg
The particular solution for an oral dose is
Eq. 28
XP,POt=KaXg(0)λ1λ20t(λ1Kpc)eλ1(tτ)+(Kpcλ2)eλ2(tτ)Kcp·eλ1(tτ)+Kcp·eλ2(tτ)eKaτdτtcp-27-43-e028.jpg
Integration of Eq. 28 becomes
Eq. 29
XP,POt=KaXg(0)λ1λ20t(λ1Kpc)eλ1te(Kaλ1)τ+(Kpcλ2)eλ2te(Kaλ2)τKcp·eλ1te(Kaλ1)τ+Kcp·eλ2te(Kaλ2)τdτtcp-27-43-e029.jpg
and
Eq. 30
XP,POt=KaXg(0)λ1λ2λ1KpcKaλ1eλ1te(Kaλ1)τ+Kpcλ2Kaλ2eλ2te(Kaλ2)τKcpKaλ1eλ1te(Kaλ1)τ+KcpKaλ2eλ2te(Kaλ2)τ0ttcp-27-43-e030.jpg
then
Eq. 31
XP,POt=KaXg(0)λ1λ2λ1KpcKaλ1eλ1teKat+Kpcλ2Kaλ2eλ2teKatKcpKaλ1eλ1teKat+KcpKaλ2eλ2teKattcp-27-43-e031.jpg
Now we can add Eq. 24, Eq. 27, and Eq. 31 for the general solution
Xt=XHt+XP,IVt+XP,POt=1λ1λ2λ1Kpceλ1t+Kpcλ2eλ2tKpc·eλ1t+Kpc·eλ2tKcp·eλ1t+Kcp·eλ2tλ1KeKcpeλ1t+Ke+Kcpλ2eλ2tX0tcp-27-43-e085.jpg
+ Rλ1λ2λ1Kpcλ1(1eλ1t)+Kpcλ2λ2(1eλ2t)Kcpλ1(1eλ1t)+Kcpλ2(1eλ2t)tcp-27-43-e086.jpg
Eq. 32
+KaXg(0)λ1λ2λ1KpcKaλ1eλ1teKat+Kpcλ2Kaλ2eλ2teKatKcpKaλ1eλ1teKat+KcpKaλ2eλ2teKattcp-27-43-e032.jpg
For intravenous only or oral administration only, replace Xg(0) or R with 0 respectively.
Special cases, when one of the divisors is 0, have different solutions.
When Kpc = Kcp = λ2 = 0 (therefore λ1 = ke), the solution reduces to that of the one-compartment model.

Three-compartment model with non-zero initial condition

The differential equation for a typical three-compartment model is
Eq. 33
dXdt=(K10+K12+K13)K21K31K12K210K130K31X+R00+KaXg(0)eKat00tcp-27-43-e033.jpg
where
K10: the elimination rate constant from the central compartment
K12: the rate constant from the central to the first peripheral compartment
K21: the rate constant from the first peripheral to the central compartment
K13: the rate constant from the central to the second peripheral compartment
K31: the rate constant from the second peripheral to the central compartment
R: infusion rate to the central compartments. Infusion rates for the peripheral compartments are 0s.
Ka: the absorption rate constant from the gut compartment
Xg(0): amount in the gut compartment at time 0. It is after the multiplication with bioavailability.
The above equation can be solved with the same way to the two-compartment model case.
Eq. 34
SIA1=1D(S+K21)(S+K31)K21(S+K31)K31(S+K21)K12(S+K31)S+K10+K12+K13S+K31K13K31K12K31K13(S+K12)K21K13S+K10+K12+K13S+K21K12K21tcp-27-43-e034.jpg
where
Eq. 35
D=S+K10+K12+K13S+K21S+K31K21K12S+K31K31K13S+K21=S3+K10+K12+K13+K21+K31S2+K10+K12+K13K21+K31+K21K31K12K21K13K31S+K10K21K31=S3+K10+K12+K13+K21+K31S2+(K10K21+K13K21+K10K31+K12K31+K21K31)S+K10K21K31tcp-27-43-e035.jpg
Let −λ1, −λ2, −λ3 be the three distinct roots of real values. (In real situations, roots are almost always three distinct real values.)
Eq. 36
SIA1=C1S+λ1+C2S+λ2+C3S+λ3C4S+λ1+C5S+λ2+C6S+λ3C7S+λ1+C8S+λ2+C9S+λ3C10S+λ1+C11S+λ2+C12S+λ3C13S+λ1+C14S+λ2+C15S+λ3C16S+λ1+C17S+λ2+C18S+λ3C19S+λ1+C20S+λ2+C21S+λ3C22S+λ1+C23S+λ2+C24S+λ3C25S+λ1+C26S+λ2+C27S+λ3tcp-27-43-e036.jpg
where
C1=(K21λ1)(K31λ1)(λ2λ1)(λ3λ1)tcp-27-43-e056.jpg C2=(K21λ2)(K31λ2)(λ1λ2)(λ3λ2)tcp-27-43-e057.jpg C3=(K21λ3)(K31λ3)(λ1λ3)(λ2λ3)tcp-27-43-e058.jpg
C4=K21(K31λ1)(λ2λ1)(λ3λ1)tcp-27-43-e059.jpg C5=K21(K31λ2)(λ1λ2)(λ3λ2)tcp-27-43-e060.jpg C6=K21(K31λ3)(λ1λ3)(λ2λ3)tcp-27-43-e061.jpg
C7=K31(K21λ1)(λ2λ1)(λ3λ1)tcp-27-43-e062.jpg C8=K31(K21λ2)(λ1λ2)(λ3λ2)tcp-27-43-e063.jpg C9=K31(K21λ3)(λ1λ3)(λ2λ3)tcp-27-43-e064.jpg
C10=K12(K31λ1)(λ2λ1)(λ3λ1)tcp-27-43-e065.jpg C11=K12(K31λ2)(λ1λ2)(λ3λ2)tcp-27-43-e066.jpg C12=K12(K31λ3)(λ1λ3)(λ2λ3)tcp-27-43-e067.jpg
C13=K1o+K12+K13λ1K31λ1K13K31(λ2λ1)(λ3λ1)tcp-27-43-e068.jpg
C14=K1o+K12+K13λ2K31λ2K13K31(λ1λ2)(λ3λ2)tcp-27-43-e069.jpg
C15=K1o+K12+K13λ3K31λ3K13K31(λ1λ3)(λ2λ3)tcp-27-43-e070.jpg
C16=K12K31(λ2λ1)(λ3λ1)tcp-27-43-e071.jpg C17=K12K31(λ1λ2)(λ3λ2)tcp-27-43-e072.jpg C18=K12K31(λ1λ3)(λ2λ3)tcp-27-43-e073.jpg
C19=K13(K21λ1)(λ2λ1)(λ3λ1)tcp-27-43-e074.jpg C20=K13(K21λ2)(λ1λ2)(λ3λ2)tcp-27-43-e075.jpg C21=K13(K21λ3)(λ1λ3)(λ2λ3)tcp-27-43-e076.jpg
C22=K21K13(λ2λ1)(λ3λ1)tcp-27-43-e077.jpg C23=K21K13(λ1λ2)(λ3λ2)tcp-27-43-e078.jpg C24=K21K13(λ1λ3)(λ2λ3)tcp-27-43-e079.jpg
C25=K1o+K12+K13λ1K21λ1K12K21(λ2λ1)(λ3λ1)tcp-27-43-e080.jpg
C26=K1o+K12+K13λ2K21λ2K12K21(λ1λ2)(λ3λ2)tcp-27-43-e081.jpg
C27=K1o+K12+K13λ3K21λ3K12K21(λ1λ3)(λ2λ3)tcp-27-43-e082.jpg
Heaviside cover-up method is convenient to get the above Cs.
Therefore, the impulse response function is
Eq. 37
expAt=C1eλ1t+C2eλ2t+C3eλ3tC4eλ1t+C5eλ2t+C6eλ3tC7eλ1t+C8eλ2t+C9eλ3tC10eλ1t+C11eλ2t+C12eλ3tC13eλ1t+C14eλ2t+C15eλ3tC16eλ1t+C17eλ2t+C18eλ3tC19eλ1t+C20eλ2t+C21eλ3tC22eλ1t+C23eλ2t+C24eλ3tC25eλ1t+C26eλ2t+C27eλ3ttcp-27-43-e037.jpg
and the homogenous solution is
Eq. 38
XHt=expAtX(0)tcp-27-43-e038.jpg
The particular solution for an intravenous infusion is
Eq. 39
XP,IVt=RC1λ1(1eλ1t)+C2λ2(1eλ2t)+C3λ3(1eλ3t)C10λ1(1eλ1t)+C11λ2(1eλ2t)+C12λ3(1eλ3t)C19λ1(1eλ1t)+C20λ2(1eλ2t)+C21λ3(1eλ3t)tcp-27-43-e039.jpg
Similarly, the particular solution for an oral dose is
Eq. 40
XP,POt=KaXg(0)C1Kaλ1(eλ1teKat)+C2Kaλ2(eλ2teKat)+C3Kaλ3(eλ3teKat)C10Kaλ1(eλ1teKat)+C11Kaλ2(eλ2teKat)+C12Kaλ3(eλ3teKat)C19Kaλ1(eλ1teKat)+C20Kaλ2(eλ2teKat)+C21Kaλ3(eλ3teKat)tcp-27-43-e040.jpg
Consequently, the summation of Eq. 38, Eq. 39, and Eq. 40 is the general solution for the three-compartment model.
Xt=XHt+XP,IVt+XP,POt=C1eλ1t+C2eλ2t+C3eλ3tC4eλ1t+C5eλ2t+C6eλ3tC7eλ1t+C8eλ2t+C9eλ3tC10eλ1t+C11eλ2t+C12eλ3tC13eλ1t+C14eλ2t+C15eλ3tC16eλ1t+C17eλ2t+C18eλ3tC19eλ1t+C20eλ2t+C21eλ3tC22eλ1t+C23eλ2t+C24eλ3tC25eλ1t+C26eλ2t+C27eλ3tX0tcp-27-43-e083.jpg
+R C1λ1(1eλ1t)+C2λ2(1eλ2t)+C3λ3(1eλ3t)C10λ1(1eλ1t)+C11λ2(1eλ2t)+C12λ3(1eλ3t)C19λ1(1eλ1t)+C20λ2(1eλ2t)+C21λ3(1eλ3t)tcp-27-43-e084.jpg
Eq. 41
+KaXg(0)C1Kaλ1(eλ1teKat)+C2Kaλ2(eλ2teKat)+C3Kaλ3(eλ3teKat)C10Kaλ1(eλ1teKat)+C11Kaλ2(eλ2teKat)+C12Kaλ3(eλ3teKat)C19Kaλ1(eλ1teKat)+C20Kaλ2(eλ2teKat)+C21Kaλ3(eλ3teKat)tcp-27-43-e041.jpg
In cases of intravenous only or oral administration only, replace Xg(0) or R with 0 respectively.
When one of the divisors becomes 0, the solution is different from the above.
When K13 = K31 = λ3 = 0, the solution reduces to that of the two-compartment model.

Conclusion

The analytical solutions of linear pharmacokinetic compartment models with a non-zero initial condition can be obtained elegantly using the matrix operations, integration, and the following facts
  • i) Pharmacokinetic compartment models for intravenous infusion or oral administration can be expressed using inhomogeneous differential equations.

  • ii) The general solution of an inhomogeneous differential equation is the sum of the homogenous solution and the particular solution.

  • iii) The particular solution is the convolution of the impulse response function from the homogenous solution and the input functions such as intravenous infusion rate or input rate from the depot compartment.

  • iv) Laplace transformation can be used to solve coupled linear differential equations.

Acknowledgments

None.

Notes

Reviewer: This article was invited and reviewed the by the editorial board members of TCP.

Conflict of interest: - Authors: Nothing to declare

- Reviewers: Nothing to declare

- Editors: Nothing to declare

TOOLS
Similar articles