2.1. Problem formulation

The transport of a solute can be modelled using the ADE, and for the contrast material used in CCTA we assume that transport is passive. This is justified when the volume fraction of the solute, or concentration, is low enough. Passive implies that the solute does not impact the material properties of the carrying fluid (here, blood), that the velocities of the solute and the carrier are equal and that no mixing of the two fluids occur. The concentration of the solute $C=C(\boldsymbol {x},t)$ is described by

(2.1)\begin{equation} \frac{\partial C}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} C = D \nabla^{2} C, \end{equation}

where $\boldsymbol {v}$ is the velocity and $D$ the diffusion constant. We introduce cylindrical coordinates $(r, \phi , z)$ with $r$ the radial, $\phi$ the circumferential and $z$ the axial coordinate and further assume that both the velocity and concentration fields are axisymmetric. The ADE (2.1) can then be written as

(2.2)\begin{equation} \frac{\partial C}{\partial t} + v_z\frac{\partial C}{\partial z} + v_r\frac{\partial C}{\partial r} = D\left(\frac{\partial^{2} C}{\partial z^{2}} + \frac{1}{r}\frac{\partial C}{\partial r} + \frac{\partial^{2} C}{\partial r^{2}}\right), \end{equation}

with the subscript indicating the corresponding component of $\boldsymbol {v}$. Finally, we assume a steady developed velocity field where $v_r = 0$ and $v_z(z, r, t) = v(r)$. This results in the following equation:

(2.3)\begin{equation} \frac{\partial C}{\partial t} + v(r)\frac{\partial C}{\partial z} = D \left(\frac{\partial^{2} C}{\partial z^{2}} + \frac{1}{r}\frac{\partial C}{\partial r} + \frac{\partial^{2} C}{\partial r^{2}}\right). \end{equation}

To make the formulation dimensionless we apply a scaling analysis using

(2.4a–e)\begin{equation} t = T \hat{t}, \quad r = a \hat{r}, \quad z = L \hat{z}, \quad v(r) = V \hat{v}(\hat{r}), \quad C(t, z, r) = C_0\hat{C}(\hat{t}, \hat{z}, \hat{r}), \end{equation}

where $T$ is the time delay between contrast arriving and reaching maximum concentration, $V$ is the cross-sectional averaged velocity, $a$ is the radius, $L$ is a characteristic length in the axial direction and $C_0$ is a characteristic concentration. With the hat to indicate a dimensionless parameter, we derive

(2.5)\begin{equation} \textit{St}\frac{\partial \hat{C}}{\partial \hat{t}} + \hat{v}(\hat{r})\frac{\partial \hat{C}}{\partial \hat{z}} = \frac{1}{\textit{Pe}}\left( \frac{1}{\hat{r}}\frac{\partial \hat{C}}{\partial \hat{r}} + \frac{\partial^{2} \hat{C}}{\partial \hat{r}^{2}} + \left(\frac{a}{L}\right)^{2} \frac{\partial^{2} \hat{C}}{\partial \hat{z}^{2}} \right), \end{equation}

where $\textit {St}$ is the Strouhal number $\textit {St} = L/V T$ and $\textit {Pe}$ the radial Pèclet number $\textit {Pe} = V a^{2}/D L$. For completeness, an estimate for each of the characteristic physiological values during CCTA is listed in table 1. For slender cylindrical domains, such as the coronary arteries we consider here, $a \ll L$, and thus we can neglect the axial diffusion in (2.5), yielding

(2.6)\begin{equation} \textit{St}\frac{\partial \hat{C}}{\partial \hat{t}} + \hat{v}(\hat{r})\frac{\partial \hat{C}}{\partial \hat{z}} = \frac{1}{\textit{Pe}}\left( \frac{1}{\hat{r}}\frac{\partial \hat{C}}{\partial \hat{r}} + \frac{\partial^{2} \mathcal{C}}{\partial \hat{r}^{2}} \right). \end{equation}

To complete the problem formulation for $\hat {C}(\hat {r},\hat {z},\hat {t})$, the equation must be supplemented by the initial and boundary conditions. Initially, i.e. for $t \leqslant 0$, no solute is present in the fluid. At the inlet, $z=0$, the input solute concentration is prescribed by the product of a time-varying function $f(t)$ and a radial variation $g(r)$. At the wall, the radial flux of $C$ is zero. This leads to the conditions

(2.7a–c)\begin{equation} \hat{C}(\hat{r},\hat{z},0)=0,\quad \hat{C}(\hat{r},0,\hat{t})=f(\hat{t})g(\hat{r}),\quad\frac{\partial \hat{C}}{\partial \hat{r}}(1,\hat{z},\hat{t})=0. \end{equation}

The axial velocity is assumed to be fully developed, implying that $v(r)$ has a parabolic profile. Choosing $V=Q/{\rm \pi} a^{2}$, the averaged velocity, with $Q$ the total flow rate, we obtain for the dimensionless velocity

(2.8)\begin{equation} \hat{v}(\hat{r}) = 2 (1- \hat{r}^{2}). \end{equation}

Unless otherwise noted, from here on we will always refer to the dimensionless formulation while omitting the hats.

Table 1.

Order estimates of the characteristic physiological values found in CCTA data.

2.2. Galerkin method

We first decompose $C(r; z; t)$ as a series of Bessel functions according to

(2.9)\begin{equation} C(r,z,t) = \sum_{m=1}^{\infty} C_m(z, t) \textrm{J}_0(\xi_m r),\quad \text{with } \textrm{J}_1(\xi_m) = 0, \end{equation}

where $\textrm {J}_0$ and $\textrm {J}_1$ are the zeroth and first-order Bessel function of the first kind with a parameter $\xi _m$ to enforce the no-flux boundary condition. Substituting this equation into (2.6), we obtain

(2.10)\begin{equation} \textit{St}\frac{\partial C_m}{\partial t}\textrm{J}_0(\xi_m r) + v(r)\frac{\partial C_m}{\partial z}\textrm{J}_0(\xi_m r) + \frac{\xi_m^{2}}{\textit{Pe}}C_m \textrm{J}_0(\xi_m r) = 0, \end{equation}

and for the boundary conditions

(2.11a,b)\begin{equation} C_m(z, 0) = 0,\quad C_m(0, t) = f(t)g_m, \end{equation}

with $g_m$ resulting from separating $g(r)$ in the same way as $C$. We employ the Galerkin method by multiplying (2.10) with $\textrm {J}_0(\xi _n r) r$ and integrating the result with respect to $r$ from $r=0$ to $1$. For practical calculations we truncate the infinite sum (2.9) at a finite number $N$, letting $m$ run from 1 to $N$. Thus, we derive

(2.12)\begin{equation} \textit{St}\frac{\partial C_n}{\partial t} + \sum_{m=1}^{N} Z_{nm}\frac{\partial C_m}{\partial z} + \frac{\xi_n^{2}}{\textit{Pe}}C_n= 0,\quad 1 \leqslant n \leqslant N, \end{equation}

with

(2.13)\begin{equation} Z_{nm} = \frac{2}{\textrm{J}_0(\xi_n)^{2}}\int_0^{1} v(r) \textrm{J}_0(\xi_m r) \textrm{J}_0(\xi_n r) r \,\textrm{d}r. \end{equation}

Substitution of $v(r) = 2(1 - r^{2})$ results for $Z_{nm}$ in

(2.14)\begin{equation} Z_{nm} = \begin{cases} 1 & \text{for } n=m=1, \\ \dfrac{4}{3} & \text{for } n=m \neq 1, \\ -8 \dfrac{(\xi_m^{2} + \xi_n^{2}) \textrm{J}_0(\xi_m)}{(\xi_m^{2} - \xi_n^{2})^{2} \textrm{J}_0(\xi_n)} & \text{for } n\neq m. \end{cases} \end{equation}

2.3. Laplace transform

We proceed by using the Laplace transform $\mathcal {L}\{C_m(z, t);t,s\} = C_{mL}(z; s)$ in time to simplify (2.12) to

(2.15a,b)\begin{equation} \sum_{m=1}^{N} Z_{nm}\frac{\partial C_{mL}}{\partial z} + \left(\frac{\xi_n^{2}}{\textit{Pe}} + \textit{St}\,s\right) C_{nL} = 0, \quad C_{nL}(0;s)=g_n f_L(s),\quad n \in (1,N), \end{equation}

with $f_L(s)$ the Laplace transform of $f(t)$. Introducing the $N$-vectors $\boldsymbol {C}_{\boldsymbol {L}} = \{ C_{nL} \}$ and $\boldsymbol {g}=\{g_n\}$ and the $N \times N$-matrices $\boldsymbol{\mathsf{Z}}=\{Z_{nm}\}$ and $\boldsymbol{\mathsf{S}}= \mathrm {Diag}\{ \xi _n^{2}/\textit {Pe} + \textit {St}\,s\}$, a diagonal matrix, we can write (2.15a,b) in matrix notation as

(2.16)\begin{equation} \boldsymbol{\mathsf{Z}}\frac{\partial \boldsymbol{C}_{L}}{\partial z} + \boldsymbol{\mathsf{S}} \boldsymbol{C}_{L} = 0. \end{equation}

The solution of the first-order ordinary differential equation (2.16) is constructed with the eigenvectors and eigenvalues of $\boldsymbol{\mathsf{Z}}^{-1} \boldsymbol{\mathsf{S}}$, resulting in

(2.17)\begin{equation} \boldsymbol{C}_L(z, s) = \sum_{k = 1}^{N} f_L(s) c_k \boldsymbol{A}_k(s) \exp({-B_k(s) z}), \end{equation}

where $\boldsymbol{A}_k$ is the $k$th eigenvector and $B_k$ is the $k$th eigenvalue of $\boldsymbol{\mathsf{Z}}^{-1} \boldsymbol{\mathsf{S}}$. Moreover, the coefficients $c_k$ are constants that can be determined by the boundary conditions

(2.18)\begin{equation} \boldsymbol{c} = \boldsymbol{\mathsf{A}}^{{-}1}\boldsymbol{g}, \end{equation}

where $\boldsymbol{\mathsf{A}}$ is the $N \times N$-matrix with columns $\boldsymbol {A}_k$, i.e. the eigensystem representation of $\boldsymbol{\mathsf{Z}}^{-1} \boldsymbol{\mathsf{S}}$. Following Abate & Whitt (Reference Abate and Whitt2006), for the inverse Laplace transform we approximate $\boldsymbol {C}$ by the finite linear combination of the transform values

(2.19)\begin{equation} \boldsymbol{C}(z, t) \approx \sum_{k = 1}^{K} \frac{\omega_k}{t}\boldsymbol{C}_L\left(z, \frac{\alpha_k}{t}\right), \end{equation}

with $K$ the number of Laplace evaluations, while $\omega$ and $\alpha$ are defined by the fixed Euler summation method for non-negative real values from Abate & Whitt (Reference Abate and Whitt2006). This approach speeds up the inversion compared to calculating the integral with a trapezoidal rule because fewer computations of $\boldsymbol {C}_L$ are needed. Note that, as stated by Abate et al., high numerical precision is needed for the inverse Laplace transform.

2.4. ADFE method

For the ADFE method we utilise the solution described in §§ 2.1–2.3 to solve the inverse problem of recovering the flow rate $Q$ given an observed concentration field $C_d$ at time of observation $t_d$. For this, ADFE only requires cross-sectional averaged concentration data $\bar {C}_d(z, t_d)$ along the vessel. Time-dependent concentration data at the inlet of the vessel needs to be available to construct $f(t)$. A schematic example of such data is shown in figure 1. The motivation for using the cross-sectional averaged concentration is attributable to the spatial resolution of CCTA; it is typically inadequate to measure the radial concentration distributions, especially near the coronary artery walls. For straight vessels, we can calculate the semi-analytical $\bar {C}_a$ as

(2.20)\begin{equation} \bar{C}_a(z,t) = 2 \int_0^{1} C(r,z,t) r \,\textrm{d}r = C_1(z,t), \end{equation}

using the computed $f(t)$ and $t=t_d$. We define the least-squares error $E_m$ between $\bar {C}_d$ and $\bar {C}_a$ as follows:

(2.21)\begin{equation} E_m = ||\bar{C}_d(z, t_d) - \bar{C}_a(z, t_d)||_2, \end{equation}

and compute $Q$ that minimises $E_m$ via an iterative minimisation method. In this work the limited-memory Broyden–Fletcher–Goldfarb–Shanno bounded algorithm available in the SciPy ‘minimise’ function, is used for the minimisation of $E_m$ (Jones, Oliphant & Peterson Reference Jones, Oliphant and Peterson2001).

Figure 1.

Schematic example of the input concentration function $f(t)$ (a) and snapshot of the cross-sectional averaged concentration at $t=1$, i.e. $\bar {C}_a(z,1)$, (b). The blue dots indicate specific concentration values evenly spaced over time and serve to show the nonlinear relationship between $f(t)$ and $\bar {C}_a(z,1)$.

2.5. Spectral-element solver

The reference solutions $C_d$ used in this work are computed using two-dimensional meshes on which (2.5), with the boundary conditions from (2.7a–c) and velocity field from (2.8) is solved numerically by the spectral-element method with high-order Legendre polynomials as the shape functions. The integration and interpolation points are chosen to be the same and are the Gauss–Lobatto–Legendre points in the $z$-direction and Gauss–Radau–Legendre points in the $r$-direction in order to avoid singularities at $r=0$ (Bernardi et al. Reference Bernardi, Dauge, Maday and Azaïez1999; Shen, Tang & Wang Reference Shen, Tang and Wang2011). Only one high-order element is used in the radial direction to ensure connectivity. For time integration the generalised-$\alpha$ method is used from Jansen, Whiting & Hulbert (Reference Jansen, Whiting and Hulbert2000). The implementation is tested through benchmark problems to ensure proper computation of solutions.