In this paper, serial theory of age-dependent compartment systems is treated. My approach seems more flexible, allowing general network structure and relaxing differentiability assumptions by basing on Riemann-Stieltjes integration.

2. Formulation of Compartment Models with Contiguous Age

A clever decomposition of the compartmental matrix $\matrix{B}$ into

\[\matrix{B} = \matrix{B}_H - \matrix{B}_{\Delta},\]

where $\matrix{B}_H$ (hollow) is $\matrix{B}$ with zeros on the diagonal and the other (diagonal) matrix complements $\matrix{B}_H$ to $\matrix{B}$, allows a compartment-age version of the McKendrick-von Foerster equation (whereas ours is a system-age version):

\begin{equation} \label{eqn:MvF} \frac{\partial \vec{m}}{\partial t} + \frac{\partial \vec{m}}{\partial \omega} = -\matrix{B}_{\Delta}(t,\omega)\,\vec{m}(t,\omega), \end{equation}

where $\omega$ is the pool age. They do not call it McKendrick-von Foerster type of equation.

  • Solving Eq. \eqref{eqn:MvF} by the method of characteristics leads to a Volterra integral equation (system) of the second kind.
  • This system is hollow in that the $i$th equation for $m_i(t,\omega)$ does not depend on $m_i(t, \omega)$ due to the zeros on the diagonal of $\matrix{B}_H$.

3. Time-Independent Coefficients

Eq. (3.7) in our notation is

\begin{equation} \label{eqn:B_ij} B_{ij}(\omega) = -\frac{1}{g_{ij}(\omega)}\,\frac{\mathrm{d}\,g_{ij}}{\mathrm{d}\,\omega} = \frac{p_{ij}(\omega)}{1-P_{ij}(\omega)}, \end{equation}

where $p_{ij}$ is the probability density function (depending on pool age) of moving from pool $j$ to pool $i$, $P_{ij}$ is the corresponding cumulative distribution function, and

\[g_{ij}(\omega) = 1 - P_{ij}(\omega)\]

is the memory function describing the fraction of mass in the compartment that remains at least until age $\omega$. Hence, the instantaneous frequency with which mass moves from $j$ to $i$ is given by

\[\frac{p_{ij}(\omega)}{g_{ij}(\omega)} = \frac{p_{ij}(\omega)}{1-P_{ij}(\omega)},\]

also known as the hazard rate in survival analysis. Hence, see Eq. \eqref{eqn:B_ij},

\[B_{ij}(\omega) = \frac{p_{ij}(\omega)}{1-P_{ij}(\omega)} = -\frac{1}{g_{ij}(\omega)}\,\frac{\mathrm{d}\,g_{ij}(\omega)}{\mathrm{d}\,\omega)}.\]

⚠️ Here they assume that mass from $j$ can only move to $i$ (serial system restriction)! We do not have that, and here it would require to handle the $P_{ij}$s more carefully with some renormalizations.

  • $P_{ij}$ can be constructed from an underlying renewal process as the cumulative distribution function of the waiting/sojourn time, $g_{ij}$ is then the survival function
  • Eq. \ref{eqn:MvF} can be stochastically interpreted as a semi-Markov process as the future only depends on the current state and the residence time (pool age) therein.
  • In general, memory functions/survival functions $g_{ij}$ can lead to hazard rates $B_{ij}$ that are not in closed form or not integrable.
  • The cumulative distribution function of leaving pool $i$ is given (in the serial as in the general case) by
\[P_i(\omega) = 1-\left[e^{-\matrix{B}^S_\Delta(\omega)}\right]_{ii},\]

where

\[\matrix{B}^S_\Delta(\omega) = \int\limits_0^\omega \matrix{B}_\Delta(t-(\omega-\omega'), \omega')\dd{\omega'}.\]

5. Demonstrations

In the advection case, the solute remains in the compartment for a finite time $T=L/v$, where $L$ is the length of the compartment (rod) and $v$ is the velocity of the fluid, and then exits to immediately enter the next compartment. Then, with

\[\Phi(x) = \mathbb{1}_{\{y:\,y\geq0\}}(x)\]

being the Heaviside step function,

\[g(\omega) = \Phi(L/v-\omega),\qquad P(\omega)=\Phi(\omega-L/v),\qquad p(\omega)=\delta(\omega-L/v).\]

7. Discussion

  • approach for more complex (non-serial), multiply connected networks is called for
  • the non-autonomous approach seems useful for climate impacts on global cycling involving diffusion
  • catalog of of age distributions would be nice
  • inverting age-structured models for calibration is an open question
  • Laplace transform approaches might be useful

Questions and Remarks

  • They derive the Volterra integral equations of the second kind by the numerical solution of the MvF PDE \eqref{eqn:MvF}. Is it true that in my document the Volterra equations appear organically simply from the theory, not from the numerical problem of solving a PDE?
    • Yes, it comes directly from the Kolmogorov forward equation for the state transition operator.
  • What are fractional (Caputo) time derivatives?
  • We do not have the restriction to serial systems.
  • We use Riemann-Stieltjes integrals, and hence we can integrate as long as we have a proper $P_{ij}$, right?
  • Am I using Laplace transforms?