# Laplace transforms for a pharmacokinetics multi-compartmental model

I am an anaesthetist trying to write some pharmacokinetics software as a pet project. Unfortunately the maths I need is a bit too much for my rusty high school calculus, and I am out of my depth. I am working with the following system of equations:

$\frac{dC_{1}}{dt}=k_{21}C_{2}(t)-k_{12}C_{1}(t)+k_{31}C_{3}(t)-k_{13}C_{1}(t)-k_{10}C_{1}(t)+C_{inf}$

$\frac{dC_{2}}{dt}=k_{12}C_{1}(t)-k_{21}C_{2}(t)$

$\frac{dC_{3}}{dt}=k_{13}C_{1}(t)-k_{31}C_{3}(t)$

$\frac{dC_{e}}{dt}=k_{e0}C_{1}(t)-k_{e0}C_{e}(t)$

These equations model the distribution of a drug in the human body according to a multi-compartmental model where:

• $C_{1}$: concentration in compartment 1 (bloodstream).
• $C_{2}$: concentration in compartment 2 (richly perfused, mostly muscle).
• $C_{3}$: concentration in compartment 3 (poorly perfused, mostly fat).
• $C_{e}$: concentration at effector site (brain).
• $C_{inf}$: concentration of infusion (for a drug given intravenously at a constant rate, the rate divided by the volume of distribution of compartment 1).
• $k_{10}$: constant of elimination from compartment 1.
• $k_{12}$, $k_{13}$, etc.: constants of equilibration between compartments.
• $k_{e0}$: constant of equilibration to effector site (it is assumed to be the same in both directions).
• t: time.

All the above are known except for $C_{1}$, $C_{2}$, $C_{3}$ and $C_{e}$. We also know what the 4 concentrations are at t = 0.

I have done Laplace transforms of these equations (with the aid of the computer program Maxima) and solved the system for the Laplace transforms. I have not been able to calculate the inverse of the transforms, so I have used Zakian's method to obtain a numerical approximation to the inverse of the Laplace transforms.

This approach works reasonably well and provides quite accurate results, but it is computationally intensive, and therefore too slow to draw graphics on the fly. Also it bugs me to use a numerical approximation when there may be an exact solution.

This link seems to indicate that for a simpler two compartment model it is possible to find an exact solution by doing Laplace transforms, solving the system of equations, and then inverting the transforms.

Is my system of equations solvable? If so, I would be very grateful for some pointers.

For practical purposes, there is no closed form solution to to this problem. A closed form solution would require you to diagonalize the coefficient matrix in a closed form. Diagonalization of an $n \times n$ matrix corresponds to solving a polynomial equation of degree $n$. For $n=2$ there is a familiar formula, the quadratic formula. For $n=3$ and $n=4$ there are cubic and quartic formulae, but they are both quite awkward and unstable; numerical root finders like Newton's method are better in practice. For $n>4$ there may be no formula in terms of radicals at all, depending on the coefficients. This might explain the discrepancy you have seen.

However, this should be quickly solvable numerically without needing to use diagonalization or Laplace transforms by simply using a standard ODE solver. A popular one is Runge-Kutta (4,5); this is a fourth order method which uses an embedded fifth order method to choose step sizes. It is implemented in MATLAB as ode45. This is a "nonstiff" method; roughly speaking, this means it is fast for "nice" problems but can be extremely slow or even fail entirely for "bad" problems. If your problem turns out to be "bad", there are standard "stiff" solvers, such as MATLAB's ode23s.

You could write your program like this:

odefun=@(C) K*C+Cinf;

tspan=[0 tmax];

C0=[1;2;3;4];

[t,C]=ode45(odefun,tspan,C0);

Here K is the coefficient matrix for your problem, and Cinf is the vector $\begin{bmatrix} C_{inf} \\ 0 \\ 0 \\ 0 \end{bmatrix}$. For illustration, the second row of $K$ is $\begin{bmatrix} k_{12} & -k_{21} & 0 & 0 \end{bmatrix}$.

• Many thanks! I don't have MATLAB, but it The Runge-Kutta method seems to be implemented in Octave too. I will try it now. If it works, I can implement it in Javascript and use it. It looks like it might be faster and more accurate than the algorithm I'm using now. – Marc Apr 11 '15 at 23:42
• I don't think diagonalization is possible at all. The constant $k_{e0}$ appears in two spots in the 4th equation, and not at all in the others, so I cannot see how you could eliminate one of those. – Marc Apr 11 '15 at 23:49
• @Marc Diagonalization is almost certainly possible in principle; "most" matrices in several precise senses are diagonalizable. I just tried it in the case where all of the $k$s are distinct, and I get a diagonalization easily using MATLAB's eig command. If you mean a closed form diagonalization, that also exists (though for certain parameter values you may get a Jordan normal form instead of a diagonalization) but as I said it is useless for all practical purposes. Out of curiosity, can you write down the values of the $k$s? – Ian Apr 11 '15 at 23:53
• Ok. Many thanks for all these explanations. All this is a bit above my head, but I will spend some time going over it and hopefully I will make some headway. The constants are different for every drug and every pharmacokinetic model. The Marsh model for propofol, for example, is used in a significant proportion of modern anaesthetics, and has the following constants: k10=0.119; k12=0.112; k13=0.042; k21=0.055; k31=0.0033; ke0=0.26 (unit is 1/min) – Marc Apr 12 '15 at 0:05
• By the way, I saw your comment to CTNT's answer. Do you think that pursuing his suggestion may give an analytical solution? – Marc Apr 12 '15 at 0:12

This should be solvable using matrix differential equation methods: basically the idea is to write the equations in the form $$\mathbf{c}'(t) = \mathbf{K} \, \mathbf{c}(t) + \mathbf{b},$$ where $\mathbf{c}(t)$ is the vector $(C_1,C_2,C_3,C_e)$, $\mathbf{K}$ is a constant matrix with components made from $k_{ij}$, to reproduce your equations and $\mathbf{b}$ is a constant vector, which in your case is $(C_{\text{inf}},0,0,0)$.

First, solve the system $$\mathbf{c}'(t) = \mathbf{K} \, \mathbf{c}(t).$$ This has solution $\mathbf{c}(t) = e^{\mathbf{K} t} \mathbf{c}(0)$, where $e^{\mathbf{K} t}$ is the matrix exponential. This calculation is most easily done by diagonalising $\mathbf{K}$, as $\mathbf{R}^{T}\mathbf{DR}$, say, where $\mathbf{D}$ is diagonal with diagonal elements $\lambda_1,\lambda_2,\lambda_3,\lambda_4$. Then $$e^{\mathbf{K} t} = \mathbf{R}^{T} e^{\mathbf{D}t}\mathbf{R},$$ where $e^{\mathbf{D}t}$ is just the diagonal matrix with diagonal elements $e^{\lambda_i t}$.

Then you need to solve the equation with the $C_{\text{inf}}$ vector in it. Here you just have to guess a solution (try a constant vector, then a constant vector times $t$ if that doesn't work, and so on).

This is certainly solvable. It is a fourth order linear system $$\dot{x}=Ax+Bu$$ with state vector $$x=\left[\matrix{C_1 & C_2 & C_3 & C_e}\right]^T$$ input $u=C_{inf}$ and matrices $A$, $B$ given by
$$A=\left[\matrix{-k_{12}-k_{13}-k_{10} & k_{21} & k_{31} & 0\\ k_{12} & -k_{21} & 0 & 0 \\ k_{13} & 0 &-k_{31} & 0 \\ k_{e0} & 0 & 0 & -k_{e0}}\right]$$ $$B=\left[\matrix{1 & 0 & 0 & 0}\right]^T$$ The Laplace transform of the state is then $$X(s)=(s\mathbb{I}_4-A)^{-1}BU(s)+(s\mathbb{I}_4-A)^{-1}x(0)$$ with state response $$x(t)=e^{At}x(0)+\int_0^t{e^{A(t-s)}Bu(s)ds}$$ Now the difficulty comes to calculating the matrix exponential. A way to solve this (among others) is through the inverse Laplace transform i.e. $$e^{At}=\mathcal{L}^{-1}\left[(s\mathbb{I}-A)^{-1}\right]$$

• Well it's the same since $\left[\matrix{1 & 0 & 0 & 0}\right]^TC_{inf}=\mathbb{I}_4\left[\matrix{C_{inf} & 0 & 0 & 0}\right]^T$ but I think my definition is more correct in the sense that there is only one input. Why do you need to append the input with extra zeros that do provide any control? – RTJ Apr 11 '15 at 23:22
• I made a simple error, sorry. (I do think this way of writing is actually a little bit confusing; here one uppercase variable is $4 \times 4$ while the other is $4 \times 1$, while one lowercase variable is $4 \times 1$ while the other is $1 \times 1$. This is good style in code, and probably standard style in control theory, but somewhat confusing to my eye in print.) – Ian Apr 11 '15 at 23:25
• I'm sorry I don't see it. $B=\left[\matrix{1 & 0 & 0 & 0}\right]^T\in\mathbb{R}^4$ and $u=C_{inf}\in\mathbb{R}$. So their product is a vector in $\mathbb{R}^4$. – RTJ Apr 11 '15 at 23:27
• Ok, no worries. – RTJ Apr 11 '15 at 23:35