Analytical solution of linear multi-compartment models with non-zero initial condition and its implementation with R

David Z. D'Argenio; Kyun-Seop Bae

doi:10.12793/tcp.2019.27.2.43

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 t_i will be calculated using the elapsed time (delta t) from the previous time point t_i-1 and an initial condition equal to the condition (amounts of each compartment) of the previous time t_i-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

X (t) = X (0) e^{{- K}_{e} t}

where X(0) is the initial condition, the amount in the body at the time 0, and Ke the elimination rate constant.

Here e^−K_et 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

\frac{d X}{d t} = - K_{e} X

If the drug is infused intravenously with the rate R (Eq. 3), the solution is Eq. 4.

Eq. 3

\frac{d X}{d t} = - K_{e} X + R

Eq. 4

X (t) = X (0) e^{{- K}_{e} t} + \frac{R}{K_{e}} (1 - e^{{- K}_{e} t})

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

\frac{d \vec{X}}{d t} = A \vec{X} + \vec{f (t)} w i t h i n i t i a l c o n d i t i o n \vec{X_{o}}

Eq. 6

\vec{X (t)} = \vec{X_{H} (t)} + \vec{X_{P} (t)}

Eq. 7

\vec{X (t)} = \exp (A t) \vec{X_{o}} + \exp (A t) * \vec{f (t)}

Eq. 8

\vec{X (t)} = \exp (A t) \vec{X_{o}} + \int_{0}^{t} \exp [A (t - τ)] f (τ) d τ

where

A: coefficient matrix

\vec{X (t)}

: general solution of the differential equation

\frac{d \vec{X}}{d t} = A \vec{X} + \vec{f (t)}

\vec{X_{H} (t)} = \exp (A t) \vec{X_{o}}

: solution of the homogenous differential equation

\frac{d \vec{X}}{d t} = A \vec{X}

exp(At): impulse response function

\vec{f (t)}

: input function

\vec{X_{P} (t)}

: particular solution of the inhomogeneous differential equation

^*: convolution operator

If a drug is infused and administered orally, the differential equation is

Eq. 9

\frac{d X}{d t} = - K_{e} X + R + K_{a} X_{g} (0) e^{- K_{a} t}

where K_a is the absorption rate constant, and X_g(0) is the orally administered dose multiplied by its bioavailability.

Then the solution is

Eq. 10

X (t) = X_{H} (t) + X_{P, I V} (t) + X_{P, P O} (t) = \exp (- K_{e} t) X (0) + \int_{0}^{t} e x p [- K_{e} (t - τ)] R d τ + \int_{0}^{t} e x p [{- K}_{e} (t - τ)] K_{a} X_{g} (0) e^{- K_{a} τ} d τ

Eq. 11

X (t) = X (0) e^{{- K}_{e} t} + \frac{R}{K_{e}} (1 - e^{{- K}_{e} t}) + \frac{K_{a} X_{g} (0)}{K_{a} - K_{e}} (e^{{- K}_{e} t} - e^{- K_{a} t})

This is interestingly the summation of each solution for intravenous bolus, infusion, and oral dosing.

If K_a = K_e = K, the solution becomes

Eq. 12

X (t) = X (0) e^{- K t} + \frac{R}{K} (1 - e^{- K t}) + X_{g} (0) K t e^{- K t}

Two-compartment model with non-zero initial condition

The differential equation for a typical two-compartment model is

Eq. 13

\frac{d \vec{X}}{d t} = [\begin{matrix} - (K_{e} + K_{c p}) & K_{p c} \\ K_{c p} & - K_{p c} \end{matrix}] \vec{X} + [\begin{matrix} R \\ 0 \end{matrix}] + [\begin{matrix} K_{a} X_{g} (0) e^{- K_{a} t} \\ 0 \end{matrix}]

where

K_e = K₁₀: the elimination rate constant from the central compartment

K_cp = K₁₂: the rate constant from the central to the peripheral compartment

K_pc = K₂₁: 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.

K_a: the absorption rate constant from the gut compartment

X_g(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

\vec{X} = {\vec{X}}_{H} (t) + {\vec{X}}_{P, I V} (t) + {\vec{X}}_{P, P O} (t) = \exp (A t) \vec{X} (0) + \int_{0}^{t} e x p [A (t - τ)] [\begin{matrix} R \\ 0 \end{matrix}] d τ + \int_{0}^{t} e x p [A (t - τ)] [\begin{matrix} K_{a} X_{g} (0) e^{- K_{a} τ} \\ 0 \end{matrix}] d τ w h e n A = [\begin{matrix} - (K_{e} + K_{c p}) & K_{p c} \\ K_{c p} & - K_{p c} \end{matrix}]

Eq. 15

\exp (A t) = L^{- 1} [{(S I - A)}^{- 1}]

L⁻¹ indicates inverse Laplace transformation.

For the matrix within the brackets,

Eq. 16

S I - A = [\begin{matrix} S + (K_{e} + K_{c p}) & - K_{p c} \\ - K_{c p} & S + K_{p c} \end{matrix}]

is inverted to have

Eq. 17

{(S I - A)}^{- 1} = \frac{[\begin{matrix} S + K_{p c} & K_{p c} \\ K_{c p} & S + (K_{e} + K_{c p}) \end{matrix}]}{S^{2} + (K_{e} + K_{c p} + K_{p c}) S + K_{e} K_{p c}}

The denominator is a quadratic equation, and if the roots are −λ₁ and −λ₂,

Eq. 18

{(S I - A)}^{- 1} = [\begin{matrix} \frac{S + K_{p c}}{(S + λ_{1}) (S + λ_{2})} & \frac{K_{p c}}{(S + λ_{1}) (S + λ_{2})} \\ \frac{K_{c p}}{(S + λ_{1}) (S + λ_{2})} & \frac{S + (K_{e} + K_{c p})}{(S + λ_{1}) (S + λ_{2})} \end{matrix}]

where

Eq. 19

λ_{1} = \frac{(K_{e} + K_{c p} + K_{p c}) + \sqrt{{(K_{e} + K_{c p} + K_{p c})}^{2} - 4 K_{e} K_{p c}}}{2}

Eq. 20

λ_{2} = \frac{(K_{e} + K_{c p} + K_{p c}) - \sqrt{{(K_{e} + K_{c p} + K_{p c})}^{2} - 4 K_{e} K_{p c}}}{2}

Here, we present only the case that roots are two distinct real values, because this is practically true in real situation.

Eq. 21

{(S I - A)}^{- 1} = [\begin{matrix} \frac{C_{1}}{S + λ_{1}} + \frac{C_{2}}{S + λ_{2}} & \frac{C_{3}}{S + λ_{1}} + \frac{C_{4}}{S + λ_{2}} \\ \frac{C_{5}}{S + λ_{1}} + \frac{C_{6}}{S + λ_{2}} & \frac{C_{7}}{S + λ_{1}} + \frac{C_{8}}{S + λ_{2}} \end{matrix}]

where

C_{1} = \frac{λ_{1} - K_{p c}}{λ_{1} - λ_{2}}

C_{2} = \frac{K_{p c} - λ_{2}}{λ_{1} - λ_{2}}

C_{3} = \frac{- K_{p c}}{λ_{1} - λ_{2}}

C_{4} = \frac{K_{p c}}{λ_{1} - λ_{2}}

C_{5} = \frac{- K_{c p}}{λ_{1} - λ_{2}}

C_{6} = \frac{K_{c p}}{λ_{1} - λ_{2}}

C_{7} = \frac{λ_{1} - K_{e} - K_{c p}}{λ_{1} - λ_{2}}

C_{8} = \frac{K_{e} + K_{c p} - λ_{2}}{λ_{1} - λ_{2}}

Therefore,

Eq. 22

\exp (A t) = [\begin{matrix} \frac{λ_{1} - K_{p c}}{λ_{1} - λ_{2}} e^{- λ_{1} t} + \frac{K_{p c} - λ_{2}}{λ_{1} - λ_{2}} e^{- λ_{2} t} & \frac{- K_{p c}}{λ_{1} - λ_{2}} e^{- λ_{1} t} + \frac{K_{p c}}{λ_{1} - λ_{2}} e^{- λ_{2} t} \\ \frac{- K_{c p}}{λ_{1} - λ_{2}} e^{- λ_{1} t} + \frac{K_{c p}}{λ_{1} - λ_{2}} e^{- λ_{2} t} & \frac{λ_{1} - K_{e} - K_{c p}}{λ_{1} - λ_{2}} e^{- λ_{1} t} + \frac{K_{e} + K_{c p} - λ_{2}}{λ_{1} - λ_{2}} e^{- λ_{2} t} \end{matrix}]

Eq.23

\exp (A t) = \frac{1}{λ_{1} - λ_{2}} [\begin{matrix} (λ_{1} - K_{p c}) e^{- λ_{1} t} + (K_{p c} - λ_{2}) e^{- λ_{2} t} & - K_{p c} \cdot e^{- λ_{1} t} + K_{p c} \cdot e^{- λ_{2} t} \\ - K_{c p} \cdot e^{- λ_{1} t} + K_{c p} \cdot e^{- λ_{2} t} & (λ_{1} - K_{e} - K_{c p}) e^{- λ_{1} t} + (K_{e} + K_{c p} - λ_{2}) e^{- λ_{2} t} \end{matrix}]

Therefore, the homogenous solution is

Eq. 24

{\vec{X}}_{H} (t) = \exp (A t) \vec{X} (0)

Then the particular solution for intravenous dosing is

Eq. 25

{\vec{X}}_{P, I V} (t) = \frac{R}{λ_{1} - λ_{2}} \int_{0}^{t} [\begin{matrix} (λ_{1} - K_{p c}) e^{- λ_{1} (t - τ)} + (K_{p c} - λ_{2}) e^{- λ_{2} (t - τ)} \\ - K_{c p} \cdot e^{- λ_{1} (t - τ)} + K_{c p} \cdot e^{- λ_{2} (t - τ)} \end{matrix}] d τ

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

{\vec{X}}_{P, I V} (t) = \frac{R}{λ_{1} - λ_{2}} {[\begin{matrix} \frac{λ_{1} - K_{p c}}{λ_{1}} e^{- λ_{1} (t - τ)} + \frac{K_{p c} - λ_{2}}{λ_{2}} e^{- λ_{2} (t - τ)} \\ \frac{- K_{c p}}{λ_{1}} e^{- λ_{1} (t - τ)} + \frac{K_{c p}}{λ_{2}} e^{- λ_{2} (t - τ)} \end{matrix}]|}_{0}^{t}

and

Eq. 27

{\vec{X}}_{P, I V} (t) = \frac{R}{λ_{1} - λ_{2}} [\begin{matrix} \frac{λ_{1} - K_{p c}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{K_{p c} - λ_{2}}{λ_{2}} (1 - e^{- λ_{2} t}) \\ \frac{- K_{c p}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{K_{c p}}{λ_{2}} (1 - e^{- λ_{2} t}) \end{matrix}]

The particular solution for an oral dose is

Eq. 28

{\vec{X}}_{P, P O} (t) = \frac{K_{a} X_{g} (0)}{λ_{1} - λ_{2}} \int_{0}^{t} [\begin{matrix} (λ_{1} - K_{p c}) e^{- λ_{1} (t - τ)} + (K_{p c} - λ_{2}) e^{- λ_{2} (t - τ)} \\ - K_{c p} \cdot e^{- λ_{1} (t - τ)} + K_{c p} \cdot e^{- λ_{2} (t - τ)} \end{matrix}] e^{- K_{a} τ} d τ

Integration of Eq. 28 becomes

Eq. 29

{\vec{X}}_{P, P O} (t) = \frac{K_{a} X_{g} (0)}{λ_{1} - λ_{2}} \int_{0}^{t} [\begin{matrix} (λ_{1} - K_{p c}) e^{- λ_{1} t} e^{- (K_{a} - λ_{1}) τ} + (K_{p c} - λ_{2}) e^{- λ_{2} t} e^{- (K_{a} - λ_{2}) τ} \\ - K_{c p} \cdot e^{- λ_{1} t} e^{- (K_{a} - λ_{1}) τ} + K_{c p} \cdot e^{- λ_{2} t} e^{- (K_{a} - λ_{2}) τ} \end{matrix}] d τ

and

Eq. 30

{\vec{X}}_{P, P O} (t) = \frac{K_{a} X_{g} (0)}{λ_{1} - λ_{2}} {[\begin{matrix} \frac{λ_{1} - K_{p c}}{{K_{a} - λ}_{1}} e^{- λ_{1} t} e^{- (K_{a} - λ_{1}) τ} + \frac{K_{p c} - λ_{2}}{{K_{a} - λ}_{2}} e^{- λ_{2} t} e^{- (K_{a} - λ_{2}) τ} \\ \frac{- K_{c p}}{K_{a} - λ_{1}} e^{- λ_{1} t} e^{- (K_{a} - λ_{1}) τ} + \frac{K_{c p}}{{K_{a} - λ}_{2}} e^{- λ_{2} t} e^{- (K_{a} - λ_{2}) τ} \end{matrix}]|}_{0}^{t}

then

Eq. 31

{\vec{X}}_{P, P O} (t) = \frac{K_{a} X_{g} (0)}{λ_{1} - λ_{2}} [\begin{matrix} \frac{λ_{1} - K_{p c}}{{K_{a} - λ}_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{K_{p c} - λ_{2}}{{K_{a} - λ}_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) \\ \frac{- K_{c p}}{K_{a} - λ_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{K_{c p}}{{K_{a} - λ}_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) \end{matrix}]

Now we can add Eq. 24, Eq. 27, and Eq. 31 for the general solution

\vec{X} (t) = {\vec{X}}_{H} (t) + {\vec{X}}_{P, I V} (t) + {\vec{X}}_{P, P O} (t) = \frac{1}{λ_{1} - λ_{2}} [\begin{matrix} (λ_{1} - K_{p c}) e^{- λ_{1} t} + (K_{p c} - λ_{2}) e^{- λ_{2} t} & - K_{p c} \cdot e^{- λ_{1} t} + K_{p c} \cdot e^{- λ_{2} t} \\ - K_{c p} \cdot e^{- λ_{1} t} + K_{c p} \cdot e^{- λ_{2} t} & (λ_{1} - K_{e} - K_{c p}) e^{- λ_{1} t} + (K_{e} + K_{c p} - λ_{2}) e^{- λ_{2} t} \end{matrix}] \vec{X} (0)

+ \frac{R}{λ_{1} - λ_{2}} [\begin{matrix} \frac{λ_{1} - K_{p c}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{K_{p c} - λ_{2}}{λ_{2}} (1 - e^{- λ_{2} t}) \\ \frac{- K_{c p}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{K_{c p}}{λ_{2}} (1 - e^{- λ_{2} t}) \end{matrix}]

Eq. 32

+ \frac{K_{a} X_{g} (0)}{λ_{1} - λ_{2}} [\begin{matrix} \frac{λ_{1} - K_{p c}}{{K_{a} - λ}_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{K_{p c} - λ_{2}}{{K_{a} - λ}_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) \\ \frac{- K_{c p}}{K_{a} - λ_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{K_{c p}}{{K_{a} - λ}_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) \end{matrix}]

For intravenous only or oral administration only, replace X_g(0) or R with 0 respectively.

Special cases, when one of the divisors is 0, have different solutions.

When K_pc = K_cp = λ₂ = 0 (therefore λ₁ = k_e), 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

\frac{d \vec{X}}{d t} = [\begin{matrix} - (K_{10} + K_{12} + K_{13}) & K_{21} & K_{31} \\ K_{12} & - K_{21} & 0 \\ K_{13} & 0 & - K_{31} \end{matrix}] \vec{X} + [\begin{matrix} R \\ 0 \\ 0 \end{matrix}] + [\begin{matrix} K_{a} X_{g} (0) e^{- K_{a} t} \\ 0 \\ 0 \end{matrix}]

where

K₁₀: the elimination rate constant from the central compartment

K₁₂: the rate constant from the central to the first peripheral compartment

K₂₁: the rate constant from the first peripheral to the central compartment

K₁₃: the rate constant from the central to the second peripheral compartment

K₃₁: 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.

K_a: the absorption rate constant from the gut compartment

X_g(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

{(S I - A)}^{- 1} = \frac{1}{D} [\begin{matrix} (S + K_{21}) (S + K_{31}) & K_{21} (S + K_{31}) & K_{31} (S + K_{21}) \\ K_{12} (S + K_{31}) & (S + K_{10} + K_{12} + K_{13}) (S + K_{31}) - K_{13} K_{31} & K_{12} K_{31} \\ K_{13} (S + K_{12}) & K_{21} K_{13} & (S + K_{10} + K_{12} + K_{13}) (S + K_{21}) - K_{12} K_{21} \end{matrix}]

where

Eq. 35

D = (S + K_{10} + K_{12} + K_{13}) (S + K_{21}) (S + K_{31}) - K_{21} K_{12} (S + K_{31}) - K_{31} K_{13} (S + K_{21}) = S^{3} + (K_{10} + K_{12} + K_{13} + K_{21} + K_{31}) S^{2} + \{(K_{10} + K_{12} + K_{13}) (K_{21} + K_{31}) + K_{21} K_{31} - K_{12} K_{21} - K_{13} K_{31}\} S + K_{10} K_{21} K_{31} = S^{3} + (K_{10} + K_{12} + K_{13} + K_{21} + K_{31}) S^{2} + (K_{10} K_{21} + K_{13} K_{21} + K_{10} K_{31} + K_{12} K_{31} + K_{21} K_{31}) S + K_{10} K_{21} K_{31}

Let −λ₁, −λ₂, −λ₃ be the three distinct roots of real values. (In real situations, roots are almost always three distinct real values.)

Eq. 36

{(S I - A)}^{- 1} = [\begin{matrix} \frac{C_{1}}{S + λ_{1}} + \frac{C_{2}}{S + λ_{2}} + \frac{C_{3}}{S + λ_{3}} & \frac{C_{4}}{S + λ_{1}} + \frac{C_{5}}{S + λ_{2}} + \frac{C_{6}}{S + λ_{3}} & \frac{C_{7}}{S + λ_{1}} + \frac{C_{8}}{S + λ_{2}} + \frac{C_{9}}{S + λ_{3}} \\ \frac{C_{10}}{S + λ_{1}} + \frac{C_{11}}{S + λ_{2}} + \frac{C_{12}}{S + λ_{3}} & \frac{C_{13}}{S + λ_{1}} + \frac{C_{14}}{S + λ_{2}} + \frac{C_{15}}{S + λ_{3}} & \frac{C_{16}}{S + λ_{1}} + \frac{C_{17}}{S + λ_{2}} + \frac{C_{18}}{S + λ_{3}} \\ \frac{C_{19}}{S + λ_{1}} + \frac{C_{20}}{S + λ_{2}} + \frac{C_{21}}{S + λ_{3}} & \frac{C_{22}}{S + λ_{1}} + \frac{C_{23}}{S + λ_{2}} + \frac{C_{24}}{S + λ_{3}} & \frac{C_{25}}{S + λ_{1}} + \frac{C_{26}}{S + λ_{2}} + \frac{C_{27}}{S + λ_{3}} \end{matrix}]

where

C_{1} = \frac{(K_{21} - λ_{1}) (K_{31} - λ_{1})}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{2} = \frac{(K_{21} - λ_{2}) (K_{31} - λ_{2})}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{3} = \frac{(K_{21} - λ_{3}) (K_{31} - λ_{3})}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

C_{4} = \frac{K_{21} (K_{31} - λ_{1})}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{5} = \frac{K_{21} (K_{31} - λ_{2})}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{6} = \frac{K_{21} (K_{31} - λ_{3})}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

C_{7} = \frac{K_{31} (K_{21} - λ_{1})}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{8} = \frac{K_{31} (K_{21} - λ_{2})}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{9} = \frac{K_{31} (K_{21} - λ_{3})}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

C_{10} = \frac{K_{12} (K_{31} - λ_{1})}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{11} = \frac{K_{12} (K_{31} - λ_{2})}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{12} = \frac{K_{12} (K_{31} - λ_{3})}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

C_{13} = \frac{(K_{1 o} + K_{12} + K_{13} - λ_{1}) (K_{31} - λ_{1}) - K_{13} K_{31}}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{14} = \frac{(K_{1 o} + K_{12} + K_{13} - λ_{2}) (K_{31} - λ_{2}) - K_{13} K_{31}}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{15} = \frac{(K_{1 o} + K_{12} + K_{13} - λ_{3}) (K_{31} - λ_{3}) - K_{13} K_{31}}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

C_{16} = \frac{K_{12} K_{31}}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{17} = \frac{K_{12} K_{31}}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{18} = \frac{K_{12} K_{31}}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

C_{19} = \frac{K_{13} (K_{21} - λ_{1})}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{20} = \frac{K_{13} (K_{21} - λ_{2})}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{21} = \frac{K_{13} (K_{21} - λ_{3})}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

C_{22} = \frac{K_{21} K_{13}}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{23} = \frac{K_{21} K_{13}}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{24} = \frac{K_{21} K_{13}}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

C_{25} = \frac{(K_{1 o} + K_{12} + K_{13} - λ_{1}) (K_{21} - λ_{1}) - K_{12} K_{21}}{(λ_{2} - λ_{1}) (λ_{3} - λ_{1})}

C_{26} = \frac{(K_{1 o} + K_{12} + K_{13} - λ_{2}) (K_{21} - λ_{2}) - K_{12} K_{21}}{(λ_{1} - λ_{2}) (λ_{3} - λ_{2})}

C_{27} = \frac{(K_{1 o} + K_{12} + K_{13} - λ_{3}) (K_{21} - λ_{3}) - K_{12} K_{21}}{(λ_{1} - λ_{3}) (λ_{2} - λ_{3})}

Heaviside cover-up method is convenient to get the above Cs.

Therefore, the impulse response function is

Eq. 37

\exp (A t) = [\begin{matrix} C_{1} e^{- λ_{1} t} + C_{2} e^{- λ_{2} t} + C_{3} e^{- λ_{3} t} & C_{4} e^{- λ_{1} t} + C_{5} e^{- λ_{2} t} + C_{6} e^{- λ_{3} t} & C_{7} e^{- λ_{1} t} + C_{8} e^{- λ_{2} t} + C_{9} e^{- λ_{3} t} \\ C_{10} e^{- λ_{1} t} + C_{11} e^{- λ_{2} t} + C_{12} e^{- λ_{3} t} & C_{13} e^{- λ_{1} t} + C_{14} e^{- λ_{2} t} + C_{15} e^{- λ_{3} t} & C_{16} e^{- λ_{1} t} + C_{17} e^{- λ_{2} t} + C_{18} e^{- λ_{3} t} \\ C_{19} e^{- λ_{1} t} + C_{20} e^{- λ_{2} t} + C_{21} e^{- λ_{3} t} & C_{22} e^{- λ_{1} t} + C_{23} e^{- λ_{2} t} + C_{24} e^{- λ_{3} t} & C_{25} e^{- λ_{1} t} + C_{26} e^{- λ_{2} t} + C_{27} e^{- λ_{3} t} \end{matrix}]

and the homogenous solution is

Eq. 38

{\vec{X}}_{H} (t) = \exp (A t) \vec{X} (0)

The particular solution for an intravenous infusion is

Eq. 39

{\vec{X}}_{P, I V} (t) = R [\begin{matrix} \frac{C_{1}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{C_{2}}{λ_{2}} (1 - e^{- λ_{2} t}) + \frac{C_{3}}{λ_{3}} (1 - e^{- λ_{3} t}) \\ \frac{C_{10}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{C_{11}}{λ_{2}} (1 - e^{- λ_{2} t}) + \frac{C_{12}}{λ_{3}} (1 - e^{- λ_{3} t}) \\ \frac{C_{19}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{C_{20}}{λ_{2}} (1 - e^{- λ_{2} t}) + \frac{C_{21}}{λ_{3}} (1 - e^{- λ_{3} t}) \end{matrix}]

Similarly, the particular solution for an oral dose is

Eq. 40

{\vec{X}}_{P, P O} (t) = K_{a} X_{g} (0) [\begin{matrix} \frac{C_{1}}{K_{a} - λ_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{C_{2}}{K_{a} - λ_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) + \frac{C_{3}}{K_{a} - λ_{3}} (e^{- λ_{3} t} - e^{- K_{a} t}) \\ \frac{C_{10}}{K_{a} - λ_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{C_{11}}{K_{a} - λ_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) + \frac{C_{12}}{K_{a} - λ_{3}} (e^{- λ_{3} t} - e^{- K_{a} t}) \\ \frac{C_{19}}{K_{a} - λ_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{C_{20}}{K_{a} - λ_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) + \frac{C_{21}}{K_{a} - λ_{3}} (e^{- λ_{3} t} - e^{- K_{a} t}) \end{matrix}]

Consequently, the summation of Eq. 38, Eq. 39, and Eq. 40 is the general solution for the three-compartment model.

\vec{X} (t) = {\vec{X}}_{H} (t) + {\vec{X}}_{P, I V} (t) + {\vec{X}}_{P, P O} (t) = [\begin{matrix} C_{1} e^{- λ_{1} t} + C_{2} e^{- λ_{2} t} + C_{3} e^{- λ_{3} t} & C_{4} e^{- λ_{1} t} + C_{5} e^{- λ_{2} t} + C_{6} e^{- λ_{3} t} & C_{7} e^{- λ_{1} t} + C_{8} e^{- λ_{2} t} + C_{9} e^{- λ_{3} t} \\ C_{10} e^{- λ_{1} t} + C_{11} e^{- λ_{2} t} + C_{12} e^{- λ_{3} t} & C_{13} e^{- λ_{1} t} + C_{14} e^{- λ_{2} t} + C_{15} e^{- λ_{3} t} & C_{16} e^{- λ_{1} t} + C_{17} e^{- λ_{2} t} + C_{18} e^{- λ_{3} t} \\ C_{19} e^{- λ_{1} t} + C_{20} e^{- λ_{2} t} + C_{21} e^{- λ_{3} t} & C_{22} e^{- λ_{1} t} + C_{23} e^{- λ_{2} t} + C_{24} e^{- λ_{3} t} & C_{25} e^{- λ_{1} t} + C_{26} e^{- λ_{2} t} + C_{27} e^{- λ_{3} t} \end{matrix}] \vec{X} (0)

+ R [\begin{matrix} \frac{C_{1}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{C_{2}}{λ_{2}} (1 - e^{- λ_{2} t}) + \frac{C_{3}}{λ_{3}} (1 - e^{- λ_{3} t}) \\ \frac{C_{10}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{C_{11}}{λ_{2}} (1 - e^{- λ_{2} t}) + \frac{C_{12}}{λ_{3}} (1 - e^{- λ_{3} t}) \\ \frac{C_{19}}{λ_{1}} (1 - e^{- λ_{1} t}) + \frac{C_{20}}{λ_{2}} (1 - e^{- λ_{2} t}) + \frac{C_{21}}{λ_{3}} (1 - e^{- λ_{3} t}) \end{matrix}]

Eq. 41

+ K_{a} X_{g} (0) [\begin{matrix} \frac{C_{1}}{K_{a} - λ_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{C_{2}}{K_{a} - λ_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) + \frac{C_{3}}{K_{a} - λ_{3}} (e^{- λ_{3} t} - e^{- K_{a} t}) \\ \frac{C_{10}}{K_{a} - λ_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{C_{11}}{K_{a} - λ_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) + \frac{C_{12}}{K_{a} - λ_{3}} (e^{- λ_{3} t} - e^{- K_{a} t}) \\ \frac{C_{19}}{K_{a} - λ_{1}} (e^{- λ_{1} t} - e^{- K_{a} t}) + \frac{C_{20}}{K_{a} - λ_{2}} (e^{- λ_{2} t} - e^{- K_{a} t}) + \frac{C_{21}}{K_{a} - λ_{3}} (e^{- λ_{3} t} - e^{- K_{a} t}) \end{matrix}]

In cases of intravenous only or oral administration only, replace X_g(0) or R with 0 respectively.

When one of the divisors becomes 0, the solution is different from the above.

When K₁₃ = K₃₁ = λ₃ = 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.