Conventional pharmacokinetic (PK) modeling of perfusion imaging typically quantifies tracer transport parameters using an isolated single-voxel approach. Although these methods are highly scalable and computationally efficient, they suffer from the fundamental assumption that each voxel acts as an isolated system with a known global inlet concentration—the arterial input function or AIF. This approximation leads to significant model errors which increase with spatial resolution (Buckley 2002, Calamante 2013, Willats and Calamante 2013, Hanson et al 2018).

In theory, this bias can be removed by the use of spatiotemporal PK models (Sourbron 2014). Implementations of this approach have mainly focused on one-compartment models with transport by diffusion (Koh 2013), convection (Zhou et al 2021, Zhang et al 2023), or both (Sourbron 2015, Elkin et al 2019, Zhang et al 2022). Hybrid approaches have also been proposed, coupling a one-compartment spatiotemporal model for interstitial transport with vascular delivery modeled by a single, global AIF (Pellerin et al 2007, Fluckiger et al 2013, Sinno et al 2021, Sainz-DeMena et al 2022, Sinno et al 2022). Experience with fully spatiotemporal two-compartment systems is extremely limited (Naevdal et al 2016).

All the above implementations apply additional constraints on the reconstructed model parameters (Shalom et al 2023 , 2024a, 2024b), for instance an assumption that diffusion is constant in space (Pellerin et al 2007), that the diffusion gradient between adjacent voxels is negligible (Fluckiger et al 2013), that parameter fields have small spatial gradients (Liu et al 2021, Zhou et al 2021, Zhang et al 2022, 2023), that transport is only radial in a lesion (Sinno et al 2021, 2022), that perfusion is modeled by Darcy flow (Naevdal et al 2016), or that parameter fields are in a known relationship to each other (Naevdal et al 2016). Constraints of this type are included to reduce the computational complexity, but it is not always clear that they are physically justified, creating a risk of new biases. For instance, it has been suggested that Darcy flow (Whitaker 1986) cannot capture the complete physiology of the capillary bed without further modification (Peyrounette et al 2018).

Previous studies did not investigate whether such additional constraints, and the biases they produce, are actually necessary. A simple parameter-counting exercise suggest they might not be: a 3D spatiotemporal model is massively overdetermined, with the data points vastly outnumbering the unknowns. The aim of this study, therefore, was to investigate whether unconstrained spatiotemporal PK models are fundamentally identifiable. The study addresses the question in 1D toy models for intravascular tracer where solutions can be generated efficiently with simple optimization routines.

The question of identifiability is investigated for two nested models of increasing complexity: a one-compartment convective blood flow model, and a two-compartment perfusion model with convective blood transport (Sourbron 2014). The perfusion model was selected as it presents challenging conditions for model fitting due to the relative similarity between the two compartments (arterial and venous).

2.1. Theory2.1.1. One-compartment blood flow model

The one compartment blood flow model is defined by the transport equation

Here 0 ≤ v ≤ 1 (dimensionless) is the blood volume fraction, and in units of mL/s/cm2 is the flow of blood (ml/s) through a unit tissue area (cm2). The local tissue concentration c (mmol/ml) measures the amount of contrast agent (mmol) per mL of blood, and is not directly measurable. Rather, what is measured in MRI is the tissue concentration C (mmol/ml), or the amount of contrast agent per mL of tissue:

We will assume throughout that blood and tissue are incompressible, so that is divergence free:

This model therefore has 3 free parameters per interior voxel, with 4 scalar fields , and one degree of freedom removed by the flow incompressibility. Expressing the transport equation in terms of the tissue concentration shows this more explicitly:

Here we introduced the blood velocity , which represents 3 degrees of freedom and is not divergence-free unless the blood volume fraction v is constant.

Without further constraints, the solution for can only be determined up to a constant: if solve the transport equations, then are solutions too, for any constant α. Therefore an additional constraint is needed to pin down the values unambiguously, for instance by assuming the volume fraction in one particular reference location is known. One approach could be to ensure that the spatial resolution is sufficiently high so that some voxels can be found which lie entirely inside a venous vessel. If these voxels are at location , we can safely assume that . Alternatively, if the velocity provides sufficient information for the particular application, and volumes or blood flows are not required, the model can be solved in the tissue concentration picture directly (equation (4)).

2.1.2. Two-compartment perfusion model

The two-compartment perfusion model involves an arterial (a) and venous (v) blood compartment, with mono-directional exchange from a to v by a perfusion field . Dropping coordinates from the definitions from here on for simplicity, the system is defined by:

where superscript indicates compartmental transport coefficients, concentrations, and volume fractions. The total volume fraction in these systems is constrained as 0 ≤ va + vv ≤ 1. Since the total blood flow is incompressible, we now have:

Incompressibility of blood flow therefore requires the total inflow to equal the total outflow across both the arterial and venous compartments. For the arterial compartment this yields an outflow characterized by the perfusion field (F) which provides a venous inflow, derived in Sourbron (2014) as:

Application of the divergence theorem converts this closed surface integral and when applied to a small volume this produces a local relation for the perfusion field (F) in terms of the arterial flow ():

The model has 7 free parameters per interior voxel: 8 scalar fields , with one degree of freedom again removed via the flow incompressibility. The measurable quantity is the total tissue concentration:

The transport equations can be written in terms of the tissue concentrations Ca = va ca and Cv = vv cv by defining arterial- and venous blood velocities and and the perfusion rate constant Kva = F/va :

This representation expresses the models directly in terms of 7 unconstrained scalar fields. As for the one-compartment case, the volumes and flows are only determined up to a constant: if and solve the equations, then and are solutions too, for any constant α. As before, the solution can be pinned down by adding a constraint such as for some suitably chosen location in a large venous vessel.

2.2. Discrete one-dimensions systems

Since the delivery of nutrients to tissue is a function of blood flow rather than blood velocity, clinical utility most likely hinges on the ability to measure flow. Hence we will simulate all systems using the representation. In order to apply the spatiotemporal compartment models to a system of N voxels measured at K time points, an upwind discretisation is applied in space and first-order discretisation in time. The result is illustrated in figure 1.

Figure 1. Illustration of the discretized compartment models. (A) A one-compartment blood flow model, showing a system with a positive flow direction (left-to-right). (B) A two-compartment perfusion model with arterial influxes and venous outfluxes at either end.

Download figure:

Standard image High-resolution image 2.2.1. One-compartment blood flow model

After discretisation, the one-compartment spatiotemporal model reduces to an N-compartment temporal model Sourbron (2014):

where Δt denotes the time step, and quantities ci (t) and vi defined at the voxel center. The rate constants kij from j to i are positive and defined by:

Here the flow fi is defined at the left interface of voxel i, and ki = ki−1,i + ki+1,i . Additional free parameters to the model are the concentrations c0(t) and cN+1(t) at the left and right boundary of the system, respectively. System influxes (J+(t) and J−(t)) are defined using these boundary concentrations (c0(t) and cN+1(t)) with the corresponding rate constants. In a 1D scenario, the incompressibility of flow implies that it is constant: fi = f. Hence the 1D one-compartment model is fully defined by the 1 + N quantities (f, vi ). For numerical stability, the time step Δt must be chosen to be smaller than the smallest voxel mean transit time:

2.2.2. Two-compartment perfusion model

After discretisation, the two-compartment spatiotemporal model becomes a system of 2N temporal compartments:

Additional free parameters to the model are the arterial- and venous concentrations and at the left and right boundary of the system, respectively. System influxes ( and ) are defined using the arterial boundary concentrations ( and ) with the corresponding rate constants. In 1D systems, the incompressibility of the flow implies that the total flow is constant () and that the arterial flow at the right boundary of a voxel is that at the left boundary minus the loss by perfusion:

This implies that the arterial flow at any boundary i is fully defined by the field Fi and the arterial flow f0a at the left boundary. For given total flow f the venous flow is then also determined everywhere (). Hence in the flow picture the discrete system is fully defined by the 3N + 2 quantities . For numerical stability, the time step must be smaller than the smallest voxel mean transit time:

2.3. Parameter reconstruction

The measured data consist of a 2D tissue concentration matrix Cmeasik with one value for each voxel i and each time point k. For given values of the discrete model parameters (volume fractions, flows and boundary concentrations), a predicted concentration Cpredik is generated by iterating the discrete equations (13) or (16), (17) with a time step Δt satisfying equations (15), (19). The resulting concentrations at high temporal resolution are then downsampled to the measured temporal resolution and scaled with the volume fractions (equation (2) or (10)).

Optimal values for the model parameters are determined by minimizing the root-mean-square difference between Cmeasik and Cpredik . The initial guesses for the total volume fraction (v) and the boundary concentrations are estimated from the data. The unknown boundary concentrations are estimated from the concentrations at the voxel nearest to the boundary. The volume fraction, up to a scaling constant, is estimated from the concentration at the last time point. Assuming a steady state has been reach at this time, tissue concentrations are directly proportional to v.

The optimization is performed iteratively over time: parameters are first optimized using only data up to an initial time t0—chosen to be after the initial peak of concentration has entered the system; subsequently the next time point is added and the parameters are optimized again, using the solutions from the previous step as initial values. This process is repeated until all time points are added.

The optimization for each time step is performed by a second-order gradient descent, after normalizing the parameters to dimensionless quantities in the range [0, 1]. For parameter values at the lower or upper bounds, a first-order method is applied. For flow values close to zero, the gradient is evaluated at zero. For a given gradient, the parameters are updated using an Adams update based linesearch (Kingma and Ba 2015). In this Adams update based approach, the moving averages of gradient and squared gradient are used to guide the optimization. The resulting update is scaled to restrict parameter updates crossing zero.

The whole pipeline from forward modeling to inversion method is implemented in python.

2.4. Simulations

The parameter reconstruction was evaluated for the one-compartment blood flow model and the two-compartment perfusion model. For each model, three digital reference objects were evaluated, detailed in tables 1 and  2). Models have total spatial dimensions of 25.6 cm with Δx = 0.8 cm, evolved to a total time of 80 s. Uniqueness and sensitivity of the solution were estimated by repeating the reconstruction with different initial guesses, noise levels and levels of temporal undersampling. Reconstruction accuracy was measured for each parameter field P by the difference between reconstruction (Prec) and ground truth (Pgt) as a percentage the mean absolute parameter value:

For comparison, the mean and standard deviation of the resulting Erel(P) distributions are reported denoted by . Since volume fractions and flows are only determined up to a constant, we measure reconstruction accuracy only for the velocities and rate constants, which are not subject to this redundancy.

Table 1. Ground truth values for the one-compartment systems. All x values used are in cm. PAIF is a population AIF (Parker et al 2006) with a defined delay (d) and a scaling factor (0 ≤ sf ≤ 1).

 Ground Truth CaseParameter123 f (ml/s/cm2)10.5−0.6 v (ml/ml) J+ (mM/s) PAIF(d = 10s, sf = 1) PAIF(d = 15s, sf = 1)0 J− (mM/s)00 PAIF(d = 15s, sf = 1)

Table 2. Ground truth values for the two-compartment systems. All x values used are in cm. λa denotes the ratio of arterial volume fraction to the total volume fraction; PAIF(d, sf ) is a population-based AIF (Parker et al 2006) with a defined delay (d) and a scaling factor (0 ≤ sf ≤ 1); G(w, h) denotes a centered Gaussian with width (w) and height (h); and Q(a, b, e) denotes a quadratic starting at a passing b at system center and ending at e.

 Ground truth caseParameter123 f a 1 (ml/s/cm2)0.90.5120.3 f v 1 (ml/s/cm2)−0.5−0.512−0.6 F (ml/s/ml) G(0.5Lx , 0.0626) −6266.25π)x) +0.021 v (ml/ml) +0.3 va + vv +0.3 λa va /v va (ml/ml) λa v ∣fa ∣/Q(19, 4.9, 19) λa v vv (ml/ml) v − va ∣fv ∣/Q(7.1, 1.5, 7.1) v − va (mM/s) PAIF(10s, 0.6) PAIF(10s, 1) PAIF(10s, 1) (mM/s) PAIF(15s, 0.4) PAIF(10s, 1) PAIF(15s, 1)

Reconstructions of noiseless data with 2 s temporal resolutions were repeated for several sets of initial values. For the one-compartment cases, these were f = ±(11, 9, 7, 4, 3, 1). For the two-compartment systems the 5 initial value sets are detailed in table 3.

Table 3. Initial guesses for parameters applied over all two-compartment cases. λa denotes the ratio of arterial volume fraction to the total volume fraction.

 Guess setParameter12345 fa 1 (ml/s/cm2)1.21.20.50.80.3 fv 1 (ml/s/cm2)−1.2−0.3−1.0−0.8−1.2 F (ml/s/ml)0.0650.0550.0450.0350.025 λa 0.40.450.50.550.6

Sensitivity to temporal undersampling was tested by reconstructing the noiseless systems with data sampled at 2, 4, 6, 8, and 10 s. Sensitivity to noise was tested by repeating reconstructions on data with signal-to-noise ratio (SNR) levels of 5, 10, 15, and 20. Gaussian noise was added with a standard deviation (σ) derived from the mean concentration:

The SNR lower limit of 5 was chosen to reflect the typical lower limit used in DCE-MRI protocols (Banerji et al 2012). For each SNR level, reconstructions are run with a given set of initial values for 5 realizations to calculate 95% confidence intervals on the reconstructed parameters.

Computations are run on a single CPU (Intel(R) Xeon(R) Gold 6152 CPU2.10 GHz), with a maximum of 10 000 iterations at each time iteration, and a gradient evaluation step 1 × 10−4 for one-compartment systems and 5 × 10−5 for two-compartment systems.

3.1. One-compartment blood flow model

Results for the noise-free one-compartment systems are summarized in figure 2, showing the solutions are accurate and independent of the initial guesses. The average error (mean ± standard deviation) across all cases is 2.9 ± 4.7% and 0.4 ± 0.3% for J and u respectively.