Laplace transforms for a pharmacokinetics multi-compartmental model

Question

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.

Answer 1

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}$.

Answer 2

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).

Answer 3

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]$$