1 Introduction
Multi-compartment models arise in many fields, for example pharmacokinetics, epidemiology, engineering, physics, biomedicine, systems theory, complexity theory and the social sciences [Reference Anderson1, Reference Godfrey2, Reference Walter and Contreras12].
Consider a compartmental system consisting of compartments numbered from 1 to n. A general compartment model can be expressed in the form [Reference Anderson1]
\begin{equation}\dot q_i(t) = I_i(t) + \sum_{\substack{j = 1 \\ j \ne i}}^n f_{i,\,j} q_j(t) - \sum_{\substack{j = 0 \\ j \ne i}}^n f_{j,i} q_i(t), \quad i = 1,2,\ldots,n,\end{equation}
where
$q_i(t)$
is the quantity of material in compartment i at time t. The rate of transfer of material from compartment j to compartment i (with
$i \ne j$
) at time t is modelled by
$f_{i,\,j} q_j(t)$
, where
$f_{i,\,j} \ge 0$
is a constant called the fractional transfer coefficient. The function
$I_i = I_i(t)$
is the rate of input of material into the ith compartment from the outside, and
$f_{0,i}$
is the fractional excretion coefficient so that
$f_{0,i} q_i(t)$
is the rate of excretion of material to the outside environment from the ith compartment at time t.
For example, when
$n = 2$
, a two-compartment model has the form
\begin{equation*}\begin{split}\dot q_1(t) & = I_1(t) + f_{1,2} q_2(t) - (f_{0,1} + f_{2,1}) q_1(t), \\\dot q_2(t) & = I_2(t) + f_{2,1} q_1(t) - (f_{0,2} + f_{1,2}) q_2(t).\end{split}\end{equation*}
A block diagram for a two-compartment system is illustrated in Figure 1.
Block diagram for a two-compartment system.
For notational convenience, define
\begin{equation*}f_{i,i} = - \sum_{\substack{j = 0 \\ j \ne i}}^n f_{j,i}, \quad i = 1,2,\ldots,n;\end{equation*}
hence the total outflow from compartment i to the other compartments and the outside environment at time t is
$f_{i,i} q_i(t)$
. In matrix-vector form, (1.1) can therefore be written as
\begin{equation}\dot q(t) = F q(t) + I(t),\end{equation}
where
\begin{equation*}q(t) = \begin{pmatrix}q_1(t) \\ \\[-7pt]q_2(t) \\ \\[-7pt]\vdots \\ \\[-7pt]q_n(t)\end{pmatrix}, \quad F =\begin{pmatrix}f_{1,1} &\quad f_{1,2} &\quad \cdots &\quad f_{1,n} \\ \\[-7pt]f_{2,1} &\quad f_{2,2} &\quad \cdots &\quad f_{2,n} \\ \\[-7pt]\vdots &\quad \vdots &\quad &\quad \vdots \\ \\[-7pt] \,f_{n,1} &\quad f_{n,2} &\quad \cdots & \quad f_{n,n}\,\end{pmatrix}, \quad I(t) = \begin{pmatrix}I_1(t) \\ \\[-7pt]I_2(t) \\ \\[-7pt]\vdots \\ \\[-7pt]I_n(t)\end{pmatrix}.\end{equation*}
Equation (1.2) is to be considered with some given initial condition q(0).
A pharmacologically relevant example that will be considered in this article for illustration purposes is a two-compartment model that describes the ingestion and subsequent metabolism of a drug taken orally. The first compartment is the gastrointestinal (GI) tract and the second compartment is the bloodstream. Let
$q_1(t)$
and
$q_2(t)$
be the drug masses at time t in the first and second compartments, respectively, while
$I_1(t)$
is the drug ingestion rate at time t. In this case,
$f_{1,2} = f_{0,1} = 0$
and
$I_2(t) = 0$
for all
$t > 0$
. Then, the pharmacokinetic model is
\begin{align}\dot q_1(t) & = I_1(t) - f_{2,1} q_1(t), \nonumber\\[3pt]\dot q_2(t) & = f_{2,1} q_1(t) - f_{0,2} q_2(t).\end{align}
A direct problem is where F, I(t) for all
$t > 0$
and q(0) are given in (1.2), and we wish to determine the solution
$q = q(t)$
for all
$t > 0$
. On the other hand, an inverse problem is where q(t) and I(t) for
$t > 0$
, as well as q(0), are given but this time we want to estimate the matrix F of fractional transfer coefficients.
For the direct problem, it is well known that the solution of (1.2) is expressed as
\begin{equation}q(t) = \textrm{e}^{t F} q(0) + \int_0^t \textrm{e}^{(t - \tau) F} I(\tau) \, \textrm{d} \tau,\end{equation}
where
$\textrm{e}^{t F}$
is the matrix exponential of F. Hence, the solution of the direct problem is tantamount to determining the matrix exponential of F. For an arbitrary matrix F, this may not be straightforward to compute although an elementary technique is due to Putzer (see [Reference Putzer6, Reference Rodrigo10], for instance).
With the aim of proposing a new method to tackle the inverse problem in this article, it is instructive to first study an important class of multi-compartment systems known as catenary models. These have the form
\begin{align}\dot q_1(t) & = I_1(t) - f_{2,1} q_1(t), \nonumber\\[3pt]\dot q_i(t) & = f_{i,i - 1} q_{i - 1}(t) - f_{i + 1,i} q_{i}(t), \quad i = 2,3,\ldots,n - 1, \nonumber\\[3pt]\dot q_n(t) & = f_{n,n-1} q_{n - 1}(t) - f_{0,n} q_n(t),\end{align}
where the compartments are arranged in a chain [Reference Anderson1]. A prototypical example of a catenary model is given in (1.3), where
$n = 2$
. For notational convenience, we identify
$f_{n + 1,n} = f_{0,n}$
, so that the catenary model (1.5) simplifies to
\begin{align}\dot q_1(t) & = I_1(t) - f_{2,1} q_1(t), \nonumber\\[3pt]\dot q_i(t) & = f_{i,i - 1} q_{i - 1}(t) - f_{i + 1,i} q_{i}(t), \quad i = 2,3,\ldots,n.\end{align}
In this paper, we utilise a unified approach via Laplace transforms to (i) solve the direct problem associated with the catenary model (1.6), (ii) use the solution of (1.6) in transform space to solve the inverse problem for the catenary model (1.6) and (iii) solve the inverse problem for the general multi-compartment model (1.2).
For (ii), we assume that the quantity of material is given in only one of the compartments and yet we determine the fractional transfer coefficients associated with all of the compartments. For example, when
$n = 2$
as in (1.3),
$I_1(t)$
is known and
$q_2(t)$
(that is the drug mass in the bloodstream) can be measured but not necessarily
$q_1(t)$
(that is the drug mass in the GI tract). The goal is to estimate
$f_{21}$
and
$f_{02}$
. We exhibit a serendipitous connection with symmetric polynomials by showing that the reciprocals of the fractional transfer coefficients are the roots of some polynomial whose coefficients are elementary symmetric polynomials. However, for (iii) where we are dealing with the general multi-compartment model (1.2), we assume that the quantities of material in all of the compartments are given and show that the matrix of fractional transfer coefficients is obtained by solving linear systems whose common coefficient matrix entries are the moments of the quantities of material.
This article is organised as follows. In Section 2, we tackle the direct problem for a catenary model with the aid of the Laplace transform. Using some of the results in Section 2, we solve the inverse problem for a catenary model in Section 3. In Section 4, we propose a parameter estimation method for the inverse problem for a general multi-compartment model. Section 5 presents the results of numerical simulations for a two-compartment catenary model. We give a brief discussion in Section 6 and conclude in Section 7.
2 Solution of the direct problem for a catenary model
Here, we solve the catenary model (1.6) using the Laplace transform. As is usual, we assume that
$q_1(0) > 0$
is given and
$q_i(0) = 0$
for
$i = 2,3,\ldots,n$
. Let
\begin{equation*}\hat q_i(s) = \mathcal L\{q_i(t);s\} = \int_0^\infty \textrm{e}^{-s t} q_i(t) \, \textrm{d} t, \quad i = 1,2,\ldots,n\end{equation*}
denote the Laplace transform of
$q_i(t)$
. Taking the Laplace transform of the first equation in (1.6) gives
\begin{equation}\hat q_1(s) = \frac{1}{s + f_{2,1}} [q_1(0) + \hat I_1(s)],\end{equation}
where
$\hat I_1(s)$
is the Laplace transform of
$I_1(t)$
. The Laplace convolution property implies that
\begin{equation}q_1(t) = q_1(0) \textrm{e}^{-f_{2,1} t} + \int_0^t \textrm{e}^{-f_{2,1} (t - \tau)} I_1(\tau) \, \textrm{d} \tau.\end{equation}
Moreover, taking the Laplace transform of the second equation in (1.6) yields
\begin{equation}s \hat q_i(s) = f_{i,i - 1} \hat q_{i - 1}(s) - f_{i + 1,i} \hat q_{i}(s) \quad \text{or} \quad \hat q_i(s) = f_{i,i - 1} \frac{\hat q_{i - 1}(s)}{s + f_{i + 1,i}}, \quad i = 2,3,\ldots,n.\end{equation}
Define the auxiliary functions
\begin{equation}F_i(s) = \prod_{k = 1}^i (s + f_{k + 1,k}), \quad \varphi_i(t) = \mathcal L^{-1}\Bigl\{\frac{1}{F_i(s)};t\Bigr\}, \quad i = 1,2,\ldots,n.\end{equation}
Note that
$F_1(s) = s + f_{2,1}$
and
$\varphi_1(t) = \textrm{e}^{-f_{2,1} t}$
. Furthermore,
$F_i(s)/F_{i - 1}(s) = s + f_{i + 1,i}$
for
$i = 2,3,\ldots,n$
. We claim that
\begin{equation}\hat q_i(s) = \frac{\prod_{k = 1}^{i - 1} f_{k + 1,k}}{F_i(s)} [q_1(0) + \hat I_1(s)], \quad i = 1,2,\ldots,n\end{equation}
satisfies (2.3). The case when
$i = 1$
is clear from (2.1). Moreover,
\begin{equation*}\frac{\hat q_i(s)}{\hat q_{i - 1}(s)} = \frac{\prod_{k = 1}^{i - 1} f_{k + 1,k}}{F_i(s)} \frac{F_{i - 1}(s)}{\prod_{k = 1}^{i - 2} f_{k + 1,k}} = \frac{f_{i,i - 1}}{s + f_{i + 1,i}}, \quad i = 2,3,\ldots,n\end{equation*}
and thus (2.5) satisfies (2.3). This proves the claim.
Hence, from the Laplace convolution property, we obtain from (2.5) that
\begin{equation}q_i(t) = \Bigl(\prod_{k = 1}^{i - 1} f_{k + 1,k}\Bigr) \Biggl[q_1(0) \varphi_i(t) + \int_0^t \varphi_i(t - \tau) I_1(\tau) \, \textrm{d} \tau\Biggr], \quad i = 1,2,\ldots,n.\end{equation}
Observe that (2.2) is included in (2.6) if we set
$i = 1$
. Although (1.6) is a particular case of (1.2), and the solution of the latter is (1.4), the solution (2.6) of (1.6) obtained via the Laplace transform is more straightforward since it makes use of the special structure of the matrix F and avoids the calculation of the matrix exponential.
Remark 2.1. It should be noted that in (2.6), we have exploited the special nearest-neighbour structure of the chain (1.6), thus avoiding the calculation of the matrix exponential.
Remark 2.2. If the fractional transfer coefficients are all distinct, then it is possible to evaluate
$\varphi_i(t)$
in (2.4) explicitly. Indeed, performing a partial fraction decomposition yields
\begin{equation*}\frac{1}{F_i(s)} = \frac{1}{\prod_{k = 1}^i (s + f_{k + 1,k})} = \sum_{k = 1}^i \frac{c_k}{s + f_{k + 1,k}}.\end{equation*}
Using L’Hôpital’s Rule, we see that
\begin{equation*}c_k = \lim_{s \rightarrow -f_{k + 1,k}} \frac{s + f_{k + 1,k}}{F_i(s)} = \frac{1}{F_i^{\prime}(-f_{k + 1,k})}.\end{equation*}
Thus,
\begin{equation*}\frac{1}{F_i(s)} = \sum_{k = 1}^i \frac{1}{F_i^{\prime}(-f_{k + 1,k})} \frac{1}{s + f_{k + 1,k}}\end{equation*}
and therefore from (2.4), we have
\begin{equation*}\varphi_i(t) = \sum_{k = 1}^i \frac{1}{F_i^{\prime}(-f_{k + 1,k})} \textrm{e}^{-f_{k + 1,k} t}, \quad F_i(s) = \prod_{k = 1}^i (s + f_{k + 1,k}), \quad i = 1,2,\ldots,n.\end{equation*}
Remark 2.3. A special case of (1.6) is when the rate of input of material into the first compartment from the outside is a sum of Dirac delta functions, that is
\begin{equation*}I_1(t) = \sum_{m = 1}^M I_{1,m} \delta(t - m T),\end{equation*}
where
$T > 0$
,
$I_{1,m} \ge 0$
for
$m = 1,2,\ldots,M$
and
$\delta$
is the Dirac delta function. In the context of pharmacokinetics, T is the dosage period,
$I_{i,m}$
is the dosage rate at time m T, and M is the number of doses after the initial dose. It follows that
\begin{equation*}\int_0^t \varphi_i(t - \tau) \delta(\tau - m T) \, \textrm{d} u = \int_{-\infty}^\infty \varphi_i(t - \tau) H(\tau) H(t - \tau) \delta(\tau - m T) \, \textrm{d} \tau = \varphi_i(t - m T) H(t - m T),\end{equation*}
where H is the usual Heaviside function and we used the property that
\begin{equation*}\int_{-\infty}^\infty g(\tau) \delta(\tau - a) \, \textrm{d} \tau = g(a), \quad a \in \mathbb R.\end{equation*}
Hence, (2.6) expresses the solution of the catenary problem (1.6) as
\begin{equation}q_i(t) = \Bigl(\prod_{k = 1}^{i - 1} f_{k + 1,k}\Bigr) \Biggl[q_1(0) \varphi_i(t) + \sum_{m = 1}^M I_{1,m} \varphi_i(t - m T) H(t - m T)\Biggr], \quad i = 1,2,\ldots,n.\end{equation}
Remark 2.4. Using Remark 2.3 in the pharmacokinetic model (1.3), straightforward calculations give
\begin{equation*}F_1(s) = s + f_{2,1}, \quad F_2(s) = (s + f_{2,1}) (s + f_{3,2}) = (s + f_{2,1}) (s + f_{0,2})\end{equation*}
from (2.4), so that
$F_1^{\prime}(s) = 1$
and
$F_2^{\prime}(s) = 2 s + f_{21} + f_{02}$
. Furthermore,
\begin{equation*}\varphi_1(t) = \textrm{e}^{-f_{2,1} t}, \quad \varphi_2(t) = \frac{1}{F_2^{\prime}(-f_{2,1})} \textrm{e}^{-f_{2,1} t} + \frac{1}{F_2^{\prime}(-f_{3,2})} \textrm{e}^{-f_{3,2} t} = \frac{1}{f_{0,2} - f_{2,1}} (\textrm{e}^{-f_{2,1} t} - \textrm{e}^{-f_{0,2} t}).\end{equation*}
From (2.7), we obtain
\begin{align}q_1(t) & = q_1(0) \textrm{e}^{-f_{2,1} t} + \sum_{m = 1}^M I_{1,m} \textrm{e}^{-f_{2,1} (t - m T)} H(t - m T), \nonumber\\[3pt] q_2(t) & = \frac{q_1(0) f_{2,1}}{f_{0,2} - f_{2,1}} \left(\textrm{e}^{-f_{2,1} t} - \textrm{e}^{-f_{0,2} t}\right) + \sum_{m = 1}^M \frac{I_{1,m} f_{2,1}}{f_{0,2} - f_{2,1}} \left[\textrm{e}^{-f_{2,1} (t - m T)} - \textrm{e}^{-f_{0,2} (t - m T)}\right] H(t - m T).\end{align}
3 Solution of the inverse problem for a catenary model
Suppose that
$q_i(t)$
is known for all
$t > 0$
for some fixed
$i \in \{1,2,\ldots,n\}$
. Also, assume that
$I_1(t)$
for all
$t > 0$
and
$q_1(0)$
are given. The task here is to determine
$f_{k + 1,k}$
for all
$k = 1,2,\ldots,i$
. In other words, only the quantity of material in the ith compartment is given but we recover all of the fractional transfer coefficients in compartments 1 to i. The idea we follow here for the inverse problem is reminiscent of the integration-based approaches to parameter estimation developed in [Reference Rodrigo and Mamon7, Reference Xi, Rodrigo and Mamon11, Reference Holder and Rodrigo3, Reference Rodrigo and Mamon8, Reference Rodrigo9, Reference Li and Rodrigo4, Reference Zulkarnaen and Rodrigo13].
For a fixed
$i \in \{1,2,\ldots,n\}$
, define
\begin{equation*}G_i(s) = \frac{\hat q_i(s)}{q_1(0) + \hat I_1(s)}.\end{equation*}
Note that
$G_i(s)$
and its derivatives with respect to s can be calculated since
$q_i(t)$
and
$I_1(t)$
are assumed to be known for all
$t > 0$
. Equations (2.5) and (2.4) give
\begin{equation}\frac{G_i(s)}{\prod_{k = 1}^{i - 1} f_{k + 1,k}} = \frac{1}{\prod_{k = 1}^{i - 1} f_{k + 1,k}} \frac{\hat q_i(s)}{q_1(0) + \hat I_1(s)} = \frac{1}{F_i(s)} = \prod_{k = 1}^i (s + f_{k + 1,k})^{-1}, \quad\end{equation}
which implies that
\begin{equation}\log(G_i(s)) - \sum_{k = 1}^{i - 1} \log(f_{k + 1,k}) = -\sum_{k = 1}^i \log(s + f_{k + 1,k}).\end{equation}
Choose a convenient value for s, say
$s = 0$
, so that
\begin{equation*}G_i(0) = \frac{1}{f_{i + 1,i}} \quad \text{or} \quad f_{i + 1,i} = \frac{1}{G_i(0)}\end{equation*}
from (3.1). Taking the mth derivative of (3.2), we see that
\begin{equation*}\frac{\textrm{d}^m}{\textrm{d} s^m} \log(G_i(s)) = (-1)^m (m - 1)! \sum_{k = 1}^i (s + f_{k + 1,k})^{-m}.\end{equation*}
At
$s = 0$
, we get
\begin{equation*}\frac{(-1)^m}{(m - 1)!} \left.\frac{\textrm{d}^m}{\textrm{d} s^m} \log(G_i(s))\right\vert_{s = 0} = \sum_{k = 1}^i \Bigl(\frac{1}{f_{k + 1,k}}\Bigr)^m = \sum_{k = 1}^{i - 1} \Bigl(\frac{1}{f_{k + 1,k}}\Bigr)^m + [G_i(0)]^m, \quad m = 1,2,\ldots.\end{equation*}
Letting
\begin{equation*}a_{i,m} = \frac{(-1)^m}{(m - 1)!} \left.\frac{\textrm{d}^m}{\textrm{d} s^m} \log(G_i(s))\right\vert_{s = 0} - [G_i(0)]^m, \quad x_k = \frac{1}{f_{k + 1,k}},\end{equation*}
we obtain the system
\begin{align*}a_{i,1} & = x_1 + x_2 + \cdots + x_{i - 1} = p_1(x_1,x_2,\ldots,x_{i - 1}), \\[3pt] a_{i,2} & = x_1^2 + x_2^2 + \cdots + x_{i - 1}^2 = p_2(x_1,x_2,\ldots,x_{i - 1}), \\[3pt] & \vdots \\[3pt] a_{i,i - 1} & = x_1^{i - 1} + x_2^{i - 1} + \cdots + x_{i - 1}^{i - 1} = p_{i - 1}(x_1,x_2,\ldots,x_{i - 1}).\end{align*}
Observe that in the above system,
$a_{i,m}$
for
$m = 1,2,\ldots,i - 1$
are known and
$p_m = p_m(x_1,x_2,\ldots,x_{i - 1})$
for
$m = 1,2,\ldots,i - 1$
are power sum symmetric polynomials [Reference Macdonald5]. Then, elementary symmetric polynomials can be expressed in terms of these power sum symmetric polynomials, that is
\begin{equation*}e_m = e_m(x_1,x_2,\ldots,x_{i - 1}) = \frac{1}{m} \sum_{k = 1}^m (-1)^{k - 1} e_{m - k}(x_1,x_2,\ldots,x_{i - 1}) p_k(x_1,x_2,\ldots,x_{i - 1}).\end{equation*}
Note that
$e_0 = e_0(x_1,x_2,\ldots,x_{i - 1}) = 1$
by definition. Hence,
\begin{equation*}e_m = \frac{1}{m} \sum_{k = 1}^m (-1)^{k - 1} e_{m - k} a_{i,k}, \quad m = 1,2,\ldots,i - 1, \quad e_0 = 1\end{equation*}
are also known. For example,
\begin{equation*}e_1 = e_0 a_{i,1} = a_{i,1}, \quad e_2 = \frac{1}{2} \left(e_1 a_{i,1} - e_0 a_{i,2}\right) = \frac{1}{2} \left(a_{i,1}^2 - a_{i,2}\right).\end{equation*}
From [Reference Macdonald5], we have that
\begin{equation}\prod_{k = 1}^{i - 1} (x - x_k) = x^{i - 1} - e_1 x^{i - 2} + e_2 x^{i - 3} + \cdots + (-1)^{i - 1} e_{i - 1}.\end{equation}
Therefore, the fractional transfer coefficients are
\begin{equation}f_{k + 1,k} = \frac{1}{x_k}, \quad k = 1,2,\ldots,i - 1, \quad f_{i + 1,i} = \frac{1}{G_i(0)},\end{equation}
where
$x_k$
for
$k = 1,2,\ldots,i - 1$
are the (positive) roots of the polynomial (3.3). This completes the solution of the inverse problem for the catenary model (1.6).
Remark 3.1. In the two-compartment model (1.3) (that is
$n =2$
), suppose that the drug mass
$q_2(t)$
in the bloodstream is given for all
$t > 0$
(that is
$i = 2$
), as well as
$q_1(0)$
and the drug ingestion rate
$I_1(t)$
for all
$t > 0$
. We want to determine the fractional transfer coefficients
$f_{2,1}$
and
$f_{0,2}$
. It follows that
\begin{equation*}G_2(s) = \frac{\hat q_2(s)}{q_1(0) + \hat I_1(s)}, \quad G_2^{\prime}(s) = \frac{[q_1(0) + \hat I_1(s)] (\hat q_2)'(s) - \hat q_2(s) (\hat I_1)'(s)}{[q_1(0) + \hat I_1(s)]^2}.\end{equation*}
As
\begin{equation*}\hat q_2(s) = \int_0^\infty \textrm{e}^{-s t} q_2(t) \, \textrm{d} t, \quad \hat I_1(s) = \int_0^\infty \textrm{e}^{-s t} I_1(t) \, \textrm{d} t,\end{equation*}
we deduce that
\begin{equation*}(\hat q_2)'(s) = -\int_0^\infty t \textrm{e}^{-s t} q_2(t) \, \textrm{d} t, \quad (\hat I_1)'(s) = -\int_0^\infty t \textrm{e}^{-s t} I_1(t) \, \textrm{d} t.\end{equation*}
Therefore,
\begin{align}\hat q_2(0) &= \int_0^\infty q_2(t) \, \textrm{d} t, \quad (\hat q_2)'(0) = -\int_0^\infty t q_2(t) \, \textrm{d} t, \nonumber\\[3pt]\hat I_1(0) &= \int_0^\infty I_1(t) \, \textrm{d} t, \quad (\hat I_1)'(0) = -\int_0^\infty t I_1(t) \, \textrm{d} t\end{align}
are known quantities (essentially the zeroth and first moments of
$q_2$
and
$I_1$
) and hence so are
\begin{equation}G_2(0) = \frac{\hat q_2(0)}{q_1(0) + \hat I_1(0)}, \quad G_2^{\prime}(0) = \frac{[q_1(0) + \hat I_1(0)] (\hat q_2)'(0) - \hat q_2(0) (\hat I_1)'(0)}{[q_1(0) + \hat I_1(0)]^2}.\end{equation}
Then, (3.4) gives
\begin{equation*}f_{2,1} = \frac{1}{x_1}, \quad f_{0,2} = f_{3,2} = \frac{1}{G_2(0)},\end{equation*}
where, from (3.3),
$x_1$
is the zero of the polynomial equation
$x - e_1 = x - a_{2,1} = 0$
or
\begin{equation*}x_1 = a_{2,1} = \frac{(-1)^1}{(1 - 1)!} \left.\frac{\textrm{d}^1}{\textrm{d} s^1} \log(G_2(s))\right\vert_{s = 0} - [G_2(0)]^1 = -\frac{G_2^{\prime}(0)}{G_2(0)} - G_2(0) = -\frac{G_2^{\prime}(0) + [G_2(0)]^2}{G_2(0)}.\end{equation*}
Thus, for the pharmacokinetic model (1.3), the fractional transfer coefficients are given by
\begin{equation}f_{2,1} = -\frac{G_2(0)}{G_2^{\prime}(0) + [G_2(0)]^2}, \quad f_{0,2} = \frac{1}{G_2(0)},\end{equation}
where
$G_2(0)$
and
$G_2^{\prime}(0)$
are as in (3.6).
4 Solution of the inverse problem for the general multi-compartment model
Let us now turn to the inverse problem for the general multi-compartment system (1.2). Here, we assume that q(t) and I(t) are known for all
$t > 0$
. Furthermore, q(0) is given. Our aim is to recover the full matrix F.
If
$\hat q(s)$
and
$\hat I(s)$
denote the (vector) Laplace transforms of q(t) and I(t), respectively, then taking the Laplace transform of (1.2) gives
\begin{equation*}s \hat q(s) - q(0) = F \hat q(s) + \hat I(s).\end{equation*}
It is not difficult to show by induction on k that
\begin{equation*}k (\hat q)^{(k - 1)}(s) + s (\hat q)^{(k)}(s) = F (\hat q)^{(k)}(s) + (\hat I)^{(k)}(s), \quad k = 1,2,\ldots,\end{equation*}
which implies that
\begin{equation}F (\hat q)^{(k)}(0) = k (\hat q)^{(k - 1)}(0) - (\hat I)^{(k)}(0), \quad k = 1,2,\ldots.\end{equation}
Define the
$n \times 1$
matrices
\begin{equation*}g_j = j (\hat q)^{(j - 1)}(0) - (\hat I)^{(j)}(0), \quad j = 1,2,\ldots,n,\end{equation*}
so that
$g_{i,j}$
denotes the entry in the ith row of
$g_j$
.
For each
$j = 1,2,\ldots,n$
, we construct an
$n \times n$
linear system from (4.1) of the form
\begin{equation*}f_{i,1} (\hat q_1)^{(j)}(0) + f_{i,2} (\hat q_2)^{(j)}(0) + \cdots + f_{i,n} (\hat q_n)^{(j)}(0) = g_{i,j}, \quad i = 1,2,\ldots,n.\end{equation*}
Hence, there are n such
$n \times n$
linear systems. Now, for each
$i = 1,2,\ldots,n$
, take the ith equation of each
$n \times n$
linear system and form a new
$n \times n$
system
\begin{equation*}f_{i,1} (\hat q_1)^{(j)}(0) + f_{i,2} (\hat q_2)^{(j)}(0) + \cdots + f_{i,n} (\hat q_n)^{(j)}(0) = g_{i,j}, \quad j = 1,2,\ldots,n.\end{equation*}
Again, there are n such
$n \times n$
linear systems. Therefore, for each
$i = 1,2,\ldots,n$
, we have the rearranged linear system
\begin{equation*}\begin{pmatrix}(\hat q_1)'(0) &\quad (\hat q_2)'(0) &\quad \cdots &\quad (\hat q_n)'(0) \\ \\[-7pt] (\hat q_1)''(0) &\quad (\hat q_2)''(0) &\quad \cdots &\quad (\hat q_n)''(0) \\ \\[-7pt]\vdots &\quad \vdots &\quad &\quad \vdots \\ \\[-7pt](\hat q_1)^{(n)}(0) &\quad (\hat q_2)^{(n)}(0) &\quad \cdots &\quad (\hat q_n)^{(n)}(0)\end{pmatrix}\begin{pmatrix}f_{i,1} \\ \\[-7pt]f_{i,2} \\ \\[-7pt]\vdots \\ \\[-7pt]\,\,f_{i,n}\end{pmatrix} =\begin{pmatrix}g_{i,1} \\ \\[-7pt]g_{i,2} \\ \\[-7pt]\vdots \\ \\[-7pt] \,\, g_{i,n}\end{pmatrix}.\end{equation*}
Solving this linear system yields the ith row of F for each
$i = 1,2,\ldots,n$
. Thus, the matrix F of fractional transfer coefficients is recovered. Note that the coefficient matrix above is the same for all
$i = 1,2,\ldots,n$
.
Remark 4.1. To exemplify the above argument, suppose that
$n = 2$
. When
$j = 1$
, (4.1) generates
\begin{align}f_{1,1} (\hat q_1)'(0) + f_{1,2} (\hat q_2)'(0) & = \hat q_1(0) - (\hat I_1)'(0) = g_{1,1}, \nonumber\\[3pt] f_{2,1} (\hat q_1)'(0) + f_{2,2} (\hat q_2)'(0) & = \hat q_2(0) - (\hat I_2)'(0) = g_{2,1},\end{align}
while when
$j = 2$
, (4.1) generates
\begin{align}f_{1,1} (\hat q_1)''(0) + f_{1,2} (\hat q_2)''(0) & = 2 (\hat q_1)'(0) - (\hat I_1)''(0) = g_{1,2}, \nonumber\\[3pt] f_{2,1} (\hat q_1)''(0) + f_{2,2} (\hat q_2)''(0) & = 2 (\hat q_2)'(0) - (\hat I_2)''(0) = g_{2,2}.\end{align}
The first equations in (4.2) and (4.3) (that is fix
$i = 1$
) together give
\begin{align*}(\hat q_1)'(0) f_{1,1} + (\hat q_2)'(0) f_{1,2} & = g_{1,1}, \\[3pt] (\hat q_1)''(0) f_{1,1} + (\hat q_2)''(0) f_{1,2} & = g_{1,2},\end{align*}
whose solution is the first row of F. Similarly, the second equations in (4.2) and (4.3) (that is fix
$i = 2$
) together provide
\begin{align*}(\hat q_1)'(0) f_{2,1} + (\hat q_2)'(0) f_{2,2} & = g_{2,1}, \\[3pt] (\hat q_1)''(0) f_{2,1} + (\hat q_2)''(0) f_{2,2} & = g_{2,2},\end{align*}
whose solution is the second row of F.
5 Numerical simulations of the inverse problem for the catenary model
For definiteness, let us consider the pharmacokinetic model (1.3). Equation (3.7) expresses the fractional transfer coefficients in terms of
$G_2(0)$
and
$G_2^{\prime}(0)$
, which in turn depend on
$\hat q_2(0)$
,
$(\hat q_2)'(0)$
,
$\hat I_1(0)$
and
$(\hat I_1)'(0)$
.
Recall that
$q_1(0)$
and
$I_1(t)$
for
$t > 0$
are assumed given. In particular, this means that
$\hat I_1(0)$
and
$(\hat I_1)'(0)$
in (3.5) are computable. However, in practice, rather than
$q_2(t)$
for all
$t > 0$
, only a finite number of measurements
$q_{2,0},q_{2,1},\ldots,q_{2,n}$
corresponding to
$t_0,t_1,\ldots,t_n$
, where
$t_0 = 0$
and
$q_{2,0} = 0$
, are available. So we approximate
\begin{align*}\hat q_2(s) &= \int_0^\infty \textrm{e}^{-s t} q_2(t) \, \textrm{d} t \simeq \int_0^{t_N} \textrm{e}^{-s t} q_2(t) \, \textrm{d} t, \\[3pt] (\hat q_2)'(s) &= -\int_0^\infty t \textrm{e}^{-s t} q_2(t) \, \textrm{d} t \simeq -\int_0^{t_N} t \textrm{e}^{-s t} q_2(t) \, \textrm{d} t,\end{align*}
where we choose
$t_N \gg t_n$
. To approximate the above integrals, we need to extrapolate the values
$q_{2,n + 1},q_{2,n + 2},\ldots,q_{2,N}$
in some way (we shall refer to this set of values as the ‘tail’ of
$q_2(t)$
). We are free to choose equally spaced points
$t_{n + 1},t_{n + 2},\ldots,t_N$
, for example.
With the help of (2.1) and (2.3), the Final Value Theorem for the Laplace transform implies that
\begin{equation*}\lim_{t \rightarrow \infty} q_1(t) = \lim_{s \rightarrow 0} s \hat q_1(s) = 0, \quad \lim_{t \rightarrow \infty} q_2(t) = \lim_{s \rightarrow 0} s \hat q_2(s) = 0\end{equation*}
and the decay is of exponential order. Thus, it is reasonable to introduce the ansatz
\begin{equation*}q_2(t) = a \textrm{e}^{-b t}, \quad t \ge t_{n - 1},\end{equation*}
where
$a, b > 0$
are to be determined. We therefore require that
$q_2(t_{n - 1}) = q_{2,n - 1}$
and
$q_2(t_n) = q_{2,n}$
to ensure continuity (note that
$t_{n - 1}$
,
$t_n$
,
$q_{2,n - 1}$
and
$q_{2,n}$
are known), that is
\begin{equation*}q_{2,n - 1} = a \textrm{e}^{-b t_{n - 1}}, \quad q_{2,n} = a \textrm{e}^{-b t_n}.\end{equation*}
From this pair of equations, we can deduce a and b, namely
\begin{equation*}b = -\frac{1}{t_n - t_{n - 1}} \log\Bigl(\frac{q_{2,n}}{q_{2,n - 1}}\Bigr), \quad a = q_{2,n} \textrm{e}^{b t_n}.\end{equation*}
To extrapolate the ‘tail’ of
$q_2(t)$
, we set
\begin{equation*}q_{2,n + 1} = a \textrm{e}^{-b t_{n + 1}}, \quad q_{2,n + 2} = a \textrm{e}^{-b t_{n + 2}}, \quad \ldots, \quad q_{2,N} = a \textrm{e}^{-b t_N}.\end{equation*}
Therefore, we have the extended data set
$\{(t_j,q_{2,j}) : j = 0,1,\ldots,n,n+1,\ldots,N\}$
, as well as
\begin{equation*}\{(t_j,\textrm{e}^{-s t_j} q_{2,j}) : j = 0,1,\ldots,n,n+1,\ldots,N\}, \quad \{(t_j,-t_j \textrm{e}^{-s t_j} q_{2,j}) : j = 0,1,\ldots,n,n+1,\ldots,N\}\end{equation*}
for any
$s \ge 0$
. We use these data sets to implement a numerical quadrature method (for example composite trapezoidal rule) to estimate
\begin{equation}\hat q_2(0) \simeq \int_0^{t_N} q_2(t) \, \textrm{d} t, \quad (\hat q_2)'(0) \simeq -\int_0^{t_N} t q_2(t) \, \textrm{d} t.\end{equation}
Together with
$q_1(0)$
,
$\hat I_1(0)$
and
$(\hat I_1)'(0)$
, we can therefore estimate
$G_2(0)$
and
$G_2^{\prime}(0)$
in (3.7), yielding the desired fractional transfer coefficients
$f_{2,1}$
and
$f_{0,2}$
for the two-compartment catenary model (1.3).
Suppose that the dosage rate function is
\begin{equation}I_1(t) = \sum_{m = 1}^M I_{1,m} \delta(t - m T).\end{equation}
Since
$\mathcal L\{\delta(t - a);s\} = \textrm{e}^{-as}$
for
$a \in \mathbb R$
, we deduce that
\begin{equation*}\hat I_1(s) = \sum_{m = 1}^M I_{1,m} \textrm{e}^{-m T s}, \quad (\hat I_1)'(s) = - T\sum_{m = 1}^M m I_{1,m} \textrm{e}^{-m T s}.\end{equation*}
Hence,
\begin{equation}\hat I_1(0) = \sum_{m = 1}^M I_{1,m}, \quad (\hat I_1)'(0) = - T\sum_{m = 1}^M m I_{1,m}.\end{equation}
We first generate theoretical data from (1.3) as follows. Take
$q_1(0) = 20$
,
$q_2(0) = 0$
,
$f_{2,1} = 2$
,
$f_{0,2} = 0.8$
,
$T = 3$
and
$M = 5$
. Absorption is typically much faster than elimination; hence,
$f_{2,1} > f_{0,2}$
. For simplicity, take
$I_{1,m} = q_1(0)$
for
$m = 1,2,\ldots,M$
. This essentially assumes that the M doses given every T days, say, are all equal to the initial dose
$q_1(0)$
. With
$n = 100$
and
$\Delta t = 18/n$
, let
$t_j = j \Delta t$
for
$j = 0,1,\ldots,n$
. Therefore, we are considering (1.3) over the time interval [0,18]. Then, we use (2.8) to find
$q_{1,j} = q_1(t_j)$
and
$q_{2,j} = q_2(t_j)$
for
$j = 0,1,\ldots,n$
. Alternatively, (1.3) can be solved numerically. The result using the analytical solution (2.8) is shown in Figure 2.
Next, we ‘keep’
$q_1(0)$
,
$I_1(t)$
for
$t > 0$
and
$\{q_{2,j} : j = 0,1,\ldots,n\}$
, and ‘hide’ everything else. We perform the ‘tail’ extrapolation procedure as explained above, choosing
$t_N = t_n + 2 t_n/3 = 5 t_n/3 = 30$
for example. The result is shown in Figure 3.
To mimic measurement errors, we add a small random perturbation (for example with a normal distribution) to the extended set
$\{q_{2,j} : j = 0,1,\ldots,n,n + 1,\ldots,N\}$
. Using the data sets
\begin{equation*}\{(t_j,q_{2,j}) : j = 0,1,\ldots,n,n+1,\ldots,N\}, \quad \{(t_j,-t_j q_{2,j}) : j = 0,1,\ldots,n,n+1,\ldots,N\},\end{equation*}
we use Scilab’s inttrap command to implement the composite trapezoidal rule and estimate the integrals
$\hat q_2(0)$
and
$(\hat q_2)'(0)$
in (5.1). Equation (5.3) is used to find
$\hat I_1(0)$
and
$(\hat I_1)'(0)$
. All of these are substituted into (3.6) to estimate
$G_2(0)$
and
$G_2^{\prime}(0)$
. Finally, we use (3.7) to calculate the fractional transfer coefficients. The result is
$f_{2,1} = 1.936521$
and
$f_{0,2} = 0.800265$
(compare with the actual values
$f_{2,1} = 2$
and
$f_{0,2} = 0.8$
used to generate the original data set).
In the above test simulation, we used
$n = 100$
for the number of measurements for
$q_2(t)$
. Of course, it is more realistic to choose a relatively small value of n. Performing the numerical simulations for different values of n produced the following results:
We can observe that estimates for
$f_{0,2}$
remain stable while those for
$f_{2,1}$
start to deviate from the correct value as n decreases. This is not unexpected as any numerical quadrature method to approximate the integrals in (5.1) becomes less accurate the coarser is the partition. In practice, we can control the coarseness of the data set over
$n + 1,n + 2,\ldots,N$
(as the ‘tail’ is extrapolated from a known decaying exponential) but not over the measured region corresponding to
$0,1,\ldots,n$
.
6 Discussion
In this section, we discuss two important issues related to the inverse problem for (1.2).
The first issue is the sensitivity of the matrix F of fractional transfer coefficients with respect to small changes in the input q. Assume that I and q(0) are fixed. Since the technique proposed here is based on integration, it is natural to quantify changes in the input q by considering changes in its moments
\begin{equation*}M_k = (-1)^k \int_0^\infty t^k q(t) \, \textrm{d} t = (\hat q)^{(k)}(0), \quad k = 0,1,\ldots.\end{equation*}
Then, the sensitivity of the fractional transfer coefficients with respect to small changes in the input can be studied by taking the Jacobian matrix of F with respect to a finite number of moments. To elucidate the above idea, consider the pharmakonetic model (1.3). In this case,
$q_2(t)$
for all
$t > 0$
is assumed to be given. The transfer coefficients
$f_{0,2}$
and
$f_{2,1}$
are as in (3.7). Defining
$c = 1/[q_1(0) + \hat I_1(0)]$
in (3.6), we can write
$G_2(0) = c M_0$
and
$G_2^{\prime}(0) = c M_1 - c^2 (\hat I_1)'(0) M_0$
. Thus, (3.7) can be expressed as
\begin{equation*}f_{0,2} = \frac{1}{c M_0}, \quad f_{2,1} = \frac{M_0}{-M_1 + c (\hat I_1)'(0) M_0 - c M_0^2}.\end{equation*}
Note that the transfer coefficients only depend on the first two moments. Straightforward differentiation yields
\begin{equation*}\frac{\partial f_{0,2}}{\partial M_0} = -\frac{1}{c M_0^2}, \quad \frac{\partial f_{0,2}}{\partial M_1} = 0\end{equation*}
and
\begin{equation*}\frac{\partial f_{2,1}}{\partial M_0} = \frac{c M_0^2 - M_1}{\left[-M_1 + c (\hat I_1)'(0) M_0 - c M_0^2\right]^2}, \quad \frac{\partial f_{2,1}}{\partial M_1} = \frac{M_0}{\left[-M_1 + c (\hat I_1)'(0) M_0 - c M_0^2\right]^2}.\end{equation*}
Therefore, an ‘integration-based’ sensitivity analysis can be performed by looking at the magnitudes of these partial derivatives.
The second issue is how to handle the case when the transfer coefficients depend on t, that is
$F = F(t)$
in (1.2). Clearly, (1.2) is not anymore solvable via the matrix exponential or the Laplace transform, and a fundamental matrix of solutions needs to be found, which is not possible in general. For the inverse problem, the idea in this paper can still be adapted if we assume that the entries of F(t) have a specific functional form. For definiteness, assume that
$f_{2,1}(t) = a_0 + a_1 t$
and
$f_{0,2}(t) = b_0 + b_1 t$
in the pharmacokinetic model (1.3), where
$a_0$
,
$a_1$
,
$b_0$
and
$b_1$
are to be determined. (More generally, we can assume that they are higher degree polynomials in t and generalise the following argument.) However, this time we have to assume that
$q_1(t)$
(and not just
$I_1(t)$
and
$q_2(t)$
) for all
$t > 0$
is known. Thus, we consider the system
\begin{equation*}\begin{split}\dot q_1(t) & = I_1(t) - (a_0 + a_1 t) q_1(t), \\[3pt] \dot q_2(t) & = (a_0 + a_1 t) q_1(t) - (b_0 + b_1 t) q_2(t),\end{split}\end{equation*}
whose Laplace transform is
\begin{align}s \hat q_1(s) - q_1(0) & = \hat I_1(s) - a_0 \hat q_1(s) - a_1 \int_0^\infty t \textrm{e}^{-s t} q_1(t) \, \textrm{d} t, \nonumber\\[3pt] s \hat q_2(s) & = a_0 \hat q_1(s) + a_1 \int_0^\infty t \textrm{e}^{-s t} q_1(t) \, \textrm{d} t - b_0 \hat q_2(s) - b_1 \int_0^\infty t \textrm{e}^{-s t} q_2(t) \, \textrm{d} t.\end{align}
Setting
$s = 0$
in (6.1) yields
\begin{align}-q_1(0) & = \hat I_1(0) - a_0 \hat q_1(0) + a_1 (\hat q_1)'(0), \nonumber\\[3pt] 0 & = a_0 \hat q_1(0) - a_1 (\hat q_1)'(0) - b_0 \hat q_2(0) + b_1 (\hat q_2)'(0).\end{align}
Differentiating (6.1) with respect to s, we have
\begin{align}\hat q_1(s) + s (\hat q_1)'(s) & = (\hat I_1)'(s) - a_0 (\hat q_1)'(s) + a_1 \int_0^\infty t^2 \textrm{e}^{-s t} q_1(t) \, \textrm{d} t, \nonumber\\[3pt] \hat q_2(s) + s (\hat q_2)'(s) & = a_0 (\hat q_1)'(s) - a_1 \int_0^\infty t^2 \textrm{e}^{-s t} q_1(t) \, \textrm{d} t - b_0 (\hat q_2)'(s) + b_1 \int_0^\infty t^2 \textrm{e}^{-s t} q_2(t) \, \textrm{d} t.\end{align}
Setting
$s = 0$
in (6.3), we obtain
\begin{align}\hat q_1(0) & = (\hat I_1)'(0) - a_0 (\hat q_1)'(0) + a_1 (\hat q_1)''(0), \nonumber\\[3pt] \hat q_2(0) & = a_0 (\hat q_1)'(0) - a_1 (\hat q_1)''(0) - b_0 (\hat q_2)'(0) + b_1 (\hat q_2)''(0).\end{align}
We then combine (6.2) and (6.4) to form a linear algebraic system for
$a_0$
,
$a_1$
,
$b_0$
and
$b_1$
, which is easily solved. In fact, the first equations in (6.2) and (6.4) give
$a_0$
and
$a_1$
, which are then substituted into the respective second equations to get
$b_0$
and
$b_1$
. Note that the coefficient matrix of the linear system involves the first three moments of
$q_1$
and
$q_2$
, which can be calculated in principle.
7 Concluding remarks
In this article, we followed a Laplace transform approach to tackle both direct and inverse problems for multi-compartment models described by systems of linear first-order ordinary differential equations.
For the direct problem, the approach is especially convenient for catenary models since it avoids the calculation of the matrix exponential. The results in Section 2 are of independent pedagogical interest since they can be used as a basis for a project for undergraduate students taking a first course in ordinary differential equations and/or mathematical modelling. The solution (1.4) of the direct problem for the general multi-compartment model (1.2) also motivates the introduction of the matrix exponential in such courses.
For the inverse problem, we investigated catenary models and obtained explicit analytical expressions for the fractional transfer coefficients in terms of elementary symmetric polynomials and the moments of the given data. We assumed that the quantity of material is given in only one compartment but were able to determine the fractional transfer coefficients in the other compartments as well. We also showed how to handle the inverse problem for a general multi-compartment model following the Laplace transform idea. However, unlike in a catenary model, in a general multi-compartment model we have to assume that the quantitities of material are available in all of the compartments so as to be able to set up the correct number of consistent linear algebraic systems. This assumption may not be realisable in practice, as we indicated even for the pharmacokinetic model (1.3).
Results of numerical simulations for catenary models by benchmarking with theoretical data (with the introduction of small random perturbations in the data to simulate real data) showed excellent results when there are many data points. As the number of data points decreases, it is expected that the accuracy will suffer as the method essentially relies on numerically integrating the data set. Since the interval of integration in the Laplace transform is a half-line, a procedure for extrapolating the ‘tail’ was derived to take into account the data outside the measured region.
A similar numerical implementation can be done for the general multi-compartment model (1.2). However, the extrapolation procedure may need to be modified depending on the matrix of fractional transfer coefficients since the exponential order at infinity may not always be true. In this general case, it is not straightforward to determine the appropriate assumptions regarding the order.
Most parameter estimation methods for multi-compartment models rely on least squares techniques (see [Reference Anderson1, Reference Godfrey2, Reference Walter and Contreras12], for instance). The Laplace transform methodology proposed in this article provides an alternative method and an additional technique for the applied mathematician’s toolbox. In a heuristic sense, instead of minimising a squared error, an integration-based approach ‘averages out the potential errors’ by taking the integrals of associated functions [Reference Zulkarnaen and Rodrigo13]. Integral transforms such as the Laplace and Mellin transforms were used in the integration-based methods proposed in [Reference Rodrigo and Mamon7, Reference Rodrigo9, Reference Xi, Rodrigo and Mamon11] because these were the appropriate integral transforms for the underlying linear differential equations. For nonlinear equations such as multi-compartment models obeying Michaelis-Menten kinetics, the ideas developed in [Reference Holder and Rodrigo3, Reference Rodrigo and Mamon8, Reference Zulkarnaen and Rodrigo13] can be adapted.
