Identifiability of spatiotemporal tissue perfusion models

Objective. Standard models for perfusion quantification in DCE-MRI produce a bias by treating voxels as isolated systems. Spatiotemporal models can remove this bias, but it is unknown whether they are fundamentally identifiable. The aim of this study is to investigate this question in silico using one-dimensional toy systems with a one-compartment blood flow model and a two-compartment perfusion model. Approach. For each of the two models, identifiability is explored theoretically and in-silico for three systems. Concentrations over space and time are simulated by forward propagation. Different levels of noise and temporal undersampling are added to investigate sensitivity to measurement error. Model parameters are fitted using a standard gradient descent algorithm, applied iteratively with a stepwise increasing time window. Model fitting is repeated with different initial values to probe uniqueness of the solution. Reconstruction accuracy is quantified for each parameter by comparison to the ground truth. Main results. Theoretical analysis shows that flows and volume fractions are only identifiable up to a constant, and that this degeneracy can be removed by proper choice of parameters. Simulations show that in all cases, the tissue concentrations can be reconstructed accurately. The one-compartment model shows accurate reconstruction of blood velocities and arterial input functions, independent of the initial values and robust to measurement error. The two-compartmental perfusion model was not fully identifiable, showing good reconstruction of arterial velocities and input functions, but multiple valid solutions for the perfusion parameters and venous velocities, and a strong sensitivity to measurement error in these parameters. Significance. These results support the use of one-compartment spatiotemporal flow models, but two-compartment perfusion models were not sufficiently identifiable. Future studies should investigate whether this degeneracy is resolved in more realistic 2D and 3D systems, by adding physically justified constraints, or by optimizing experimental parameters such as injection duration or temporal resolution.


Introduction
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(Shalom et al , 2024a(Shalom et al , 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(Sinno et al , 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.

Methodology
The question of identifiability is investigated for two nested models of increasing complexity: a onecompartment 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).

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 f  in units of mL/s/cm 2 is the flow of blood (ml/s) through a unit tissue area (cm 2 ).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 f  is divergence free: This model therefore has 3 free parameters per interior voxel, with 4 scalar fields v f , ( )


, 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 u f v   = , which represents 3 degrees of freedom and is not divergencefree unless the blood volume fraction v is constant.
Without further constraints, the solution for v f , ( ) )  a a 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 x 0  , we can safely assume that v x 1 0 ( )  = .Alternatively, if the velocity u  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)).

Two-compartment perfusion model
The two-compartment perfusion model involves an arterial (a) and venous (v) blood compartment, with monodirectional exchange from a to v by a perfusion field F x ( )  .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 v a + v v 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 (f a  ): The model has 7 free parameters per interior voxel: 8 scalar fields v f v f , ; , , 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 C a = v a c a and C v = v v c v by defining arterial-and venous blood velocities u f v a a a and the perfusion rate constant This representation expresses the models directly in terms of 7 unconstrained scalar fields.As for the onecompartment case, the volumes and flows are only determined up to a constant: if v f , a a ( )  a a are solutions too, for any constant α.As before, the solution can be pinned down by adding a constraint such as v x 1 v 0 ( )  = for some suitably chosen location x 0  in a large venous vessel.

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 v f , ( )  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.

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 c i (t) and v i defined at the voxel center.The rate constants k ij from j to i are positive and defined by: Here the flow f i is defined at the left interface of voxel i, and k i = k i−1,i + k i+1,i .Additional free parameters to the model are the concentrations c 0 (t) and c N+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 (c 0 (t) and c N+1 (t)) with the corresponding rate constants.In a 1D scenario, the incompressibility of flow implies that it is constant: Hence the 1D one-compartment model is fully defined by the 1 + N quantities ( f, v i ).For numerical stability, the time step Δt must be chosen to be smaller than the smallest voxel mean transit time:

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 c t c t , at the left and right boundary of the system, respectively.System influxes (J t a ( ) ) are defined using the arterial boundary concentrations (c t a 0 ( ) and c t N a 1 ( ) ) 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 F i and the arterial flow f 0 a at the left boundary.For given total flow f the venous flow is then also determined everywhere (f Hence in the flow picture the discrete system is fully defined by the 3N For numerical stability, the time step must be smaller than the smallest voxel mean transit time:

Parameter reconstruction
The measured data consist of a 2D tissue concentration matrix C meas ik 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 C pred ik 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 C meas ik and C pred ik .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 t 0 -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.

Simulations
The parameter reconstruction was evaluated for the one-compartment blood flow model and the twocompartment 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 (P rec ) and ground truth (P gt ) as a percentage the mean absolute parameter value:   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.
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.

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 E rel ¯(mean ± standard deviation) across all cases is 2.9 ± 4.7% and 0.4 ± 0.3% for J and u respectively.
Concentration-time data reconstructed from recovered parameter values at one of the initial guesses are shown in figure 3. Deviations between the ground truth concentration and recovered parameter profiles are not visually detectable and are between ±2.5% of the maximal concentration in each case.
The effect of SNR and sampling interval on parameter reconstruction is summarized in figure 4. Subjectlevel results are included in the supplementary information (figures S1 and S2).Average E rel ¯across all simulations at SNR 5 is 2.4 ± 2.8% and 16.1 ± 15.6% for u and J, respectively.The results show the expected behavior with increasing accuracy and precision at higher SNR and smaller Dt in all parameters.Velocities are substantially more robust to noise than the influxes, and more accurate and precise at smaller Dt levels.

Two-compartment perfusion model
Results for the noise-free two-compartment systems are summarized in figure 5, showing that the solutions are generally less well determined than in the one-compartment case.Reconstruction of arterial velocity (u a ) and influx (J a ) is most accurate with lowest E rel ¯of 27.1 ± 82.0% and 13.0 ± 38.1%, respectively.The perfusion rate (K va ) and venous velocity (u v ) are least precise with E rel ¯of 44.7% ± 76.5%, and 54.9% ± 121.4% across all cases, respectively.
Ground-truth and reconstructed concentrations for initial guess 4 and noise-free data are shown in figure 6.Despite the substantial errors in the reconstructed parameters, the reconstructed concentration is close to the ground truth and visually virtually indistinguishable.Since an accurate fit to the data is obtained with inaccurate parameters, this shows that multiple solutions are compatible with the observations.Figure 7(a) shows the impact of SNR and undersampling on parameter accuracy and precision, showing the expected trend of increasing accuracy and precision at higher SNR and smaller Dt.The magnitude of the error is generally comparable between parameters except for the perfusion rate K va , which is more sensitive to noise than the other parameters, and the venous velocity u v , which appears particularly sensitive to undersampling.Supporting figures S3, S4, S5, S6, S7, S8 illustrate these effects in more detail for all three example cases.

Discussion
The aim of this in-silico study was to determine if unconstrained spatiotemporal models for DCE-MRI are fundamentally identifiable.The data indicate that this is the case for one-compartmental blood flow models, but not for two-compartmental perfusion models.
In the absence of significant measurement error, parameters of the one-compartment model can be reconstructed accurately without imposing additional constraints on the model.They are identical even with widely different choices of the initial guesses, suggesting the solution is also unique.Reconstructions of the influxes at the boundary of the system are also accurate, confirming the idea that spatiotemporal models remove the need for a separate measurement of an arterial input function.
While the analysis in this study used non-linear optimization, the uniqueness of the one-compartmental solutions aligns with the fact that the model equations can be recast in a linear form (Sourbron 2014)-in a similar way as for standard temporal one-compartmental models (Flouri et al 2016).We chose not to implement the model in the linear form as this is known to be more noise-sensitive, and does not translate as easily to the two-compartmental scenario.Measurement error (noise and undersampling) naturally reduces the accuracy and precision of the parameters, but in a predictable and expected manner.One open question is how the parameter accuracy and precision compares to a conventional analysis using measured input functions and standard voxel-by-voxel temporal models.
Parameter reconstructions are significantly less accurate for the two-compartmental perfusion model.Multiple solutions have been found that are compatible with the data, and therefore a single optimal solution cannot be identified using a goodness-of-fit criterion alone.One approach to resolving the degeneracy in the two-compartmental model may involve modifying the experimental conditions to increase the structure in the data.However, the options in DCE-MRI are limited.The smallest sampling interval considered in this study was 2 s, and therefore there may be some room for improvement by considering faster scan sequences.Beyond that, the only additional variable that can be modified substantially is the injection protocol.The current setup uses a single bolus injection, and while this can easily be modified in practice to split the dose over two injections (Ingrisch and Sourbron 2013), there is currently little evidence that this translates to more accurate solutions.Hence this may indicate that additional constraints are needed to pin down multi-compartmental spatiotemporal models.Possible solutions previously proposed for one-compartment systems may well translate to two-compartment systems, such as the use of Darcy flows or other physical constraints to reduce the number of free variables (Naevdal et al 2016), adding regularization to impose smoothness of the solution  (Liu et al 2021, Zhou et al 2021, Zhang et al 2022, 2023), fixing less critical parameters to literature values (Pellerin et al 2007), or reverting to a measured AIF at the boundaries of the imaging slab.The use of physical constraints derived from principles of fluid dynamics and porous media theory presents a particularly attractive approach as it also provides a mechanism for studying the mechanical properties of physiological flow.While such constraints may not be necessary for one-compartmental systems, they may prove essential in the multicompartment case.
Beyond modifying experimental conditions or imposing additional constraints, another strategy for reducing the degeneracy in the solutions may well be to improve the optimization itself.Setting suitable initial conditions, for instance, may well help to bias the solution towards the correct value, and may be feasible without loosing generality.For instance, exploratory simulations with the 1D toy models suggest that initial values where the arterial velocity is higher than the venous velocity leads to better parameter recovery than randomly chosen initial values, and this is consistent with physical reality.An alternative approach, common for instance in other inverse problems in imaging such as coregistration (Studholme et al 1996, Maes et al 1999), may be to employ a multi-resolution approach, fitting parameters initially at coarse resolution and then stepwise refining the estimates until the image resolution is reached.Additionally, considering the observation that estimates are most accurate in the arterial parameters, an improvement may be possible by reparametrizing the model in terms of the arterial flow field f a (x) rather than using the perfusion field F(x) as a primary variable.Furthermore, volume fractions based on spatial variation reported in vivo could be applied.Although, an investigation into the spatial variation of volume fractions showed that deviation from the ground truth concentration and volume fractions varied similarly in areas of high and low spatial volume fraction variation.Finally, solutions proposed for temporal model fitting in DCE-MRI may well help in spatiotemporal modeling as well, such as the use of model selection, which can potentially be generalized to a voxel-by-voxel approach, and/or using the results of one-compartment fits to initialize a two-compartment analysis.The optimal strategy may also depend on the parameter that is the primary interest of the measurement.As shown, results are considerably more reliable in the properties of the upstream (arterial) compartment compared to the distal (venous) compartment and particularly the exchange parameter (perfusion) itself.The venous compartment is downstream and determined by perfusion from the upstream artery.Therefore, errors may be compounded, or the optimization hindered due to the interplay of these many parameters, potentially increasing the sensitivity to measurement error.Hence in clinical applications where the primary aim is to characterize the arterial system, issues of uniqueness identified in this study may be less critical.Unfortunately, the interest in many key clinical applications of perfusion imaging, such as acute stroke (Demeestere et al 2020) or cancer (van Dijken et al 2019), is primarily in a measurement of perfusion as this is a key metric to understanding tissue viability or metabolic activity.
This article considers purely intra-vascular models with either a singular vascular compartment or separate arterial and venous compartments with direct exchange from perfusion (equation ( 9)).This is in contrast to the widely applied single voxel approaches of the tofts-kety model (TKM) and extended tofts-kety model (ETKM) which describe tissue in terms of a vascular and an extravascular compartment.To probe perfusion the blood flow should be considered which extends the ETKM to the standard two-compartment exchange model (Sourbron and Buckley 2012).To assess permeability, the TKM and ETKM mainly utilize the K trans parameter which characterizes the transfer of contrast agent into the extravascular space from the vascular space.An over estimation of parameters from the TKM has been reported (Sinno et al 2021) due to the effect of inter-voxel exchange which is neglected.While the result in this study has shown that measurement of perfusion comes with significant numerical error, the use of spatiotemporal modeling does remove the equally substantial error that comes from assuming a single upstream feeding artery (Calamante et al 2006).It is currently unknown whether, and to what extent, this offsets the numerical reconstruction errors observed in the spatiotemporal model.
This study is obviously limited by the use a of a one-dimensional toy model.In reality, fully unconstrained spatiotemporal modeling for DCE is only relevant when applied to 3D data, as through-plane exchange of indicator cannot be excluded in realistic scenarios.However, application of multi-compartmental modeling in 3D comes with significant computational challenges that are currently largely unresolved.Standard gradientdescent type optimization as performed in this study is unlikely to be practically feasible in 3D, though this has not yet been fully explored.The use of 1D models allows for a flexible exploration of fundamental issues of parameter identifiability, but there is no guarantee that the findings translate to the 3D scenario.Indeed, 3D data are significantly more entangled due to the spatial connections in the other dimensions, and this may well help to resolved any degeneracies found in 1D.Future studies should therefore focus in the first place on developing computational methods that are able to solve spatiotemporal two-compartment models in reasonable computation times, before the issue of parameter identifiability can be investigated in-silico in 3D data.Recent developments in deep learning, specifically the use of physics informed neural networks (PINNS) and their successful application in related problems, has offered some hope that a solution may be technically feasible.

Conclusions
This study provides proof of concept that one-compartmental blood flow models are fully identifiable and do not require a separate measurement of the AIF.Arterial properties of two-compartmental perfusion models have comparable accuracy but perfusion fields and venous flows cannot be measured reliably.Future studies should focus on exploring the use of physical constraints, improved optimization and on development of computational solutions for the 3D case.

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.

Figure 2 .
Figure 2. Parameter reconstructions for all noise-free one-compartment system cases.The solid black line indicates the ground truth, and the colored dots show reconstructions with different initial guesses.

Figure 3 .
Figure 3.Comparison of the recovered concentration values from the retrieved parameters against the ground truth, alongside the percentage difference for the maximal concentration.Shown are one-compartment cases 1, 2 and 3 in rows (a), (b), and (c), respectively.Differences above or below ±2.5% of the maximum concentration value are shown by dark red or dark blue, respectively.

Figure 4 .
Figure 4. Box plots of reconstruction errors for all one-compartment parameters as a function of SNR (a) and temporal sampling Dt (b).

Figure 5 .
Figure 5. Parameter reconstructions for all noise-free two-compartment system cases.The solid black line indicates the ground truth, and the colored dots show reconstructions with different initial guesses.

Figure 6 .
Figure 6.Comparison of the recovered concentration values from the retrieved parameters against the ground truth, alongside the percentage difference for the maximal concentration.Shown are two-compartment cases 1, 2 and 3 in rows (a), (b), and (c), respectively for initial guess 4. Differences above or below ±5% of the maximum concentration value are shown by dark red or dark blue, respectively.

Figure 7 .
Figure7.Distribution of each parameter error for the two-compartment system relative to the absolute mean parameter value within each system.The distribution is shown across all 3 cases for (a) all noise realizations at each SNR; and (b) undersampling rates.

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

Table 2 .
(Parker et al 2006) 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; P AIF (d, s f ) is a population-based AIF(Parker et al 2006)with a defined delay (d) and a scaling factor (0 s f 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.
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 3 .
Initial guesses for parameters applied over all two-compartment cases.λ a denotes the ratio of arterial volume fraction to the total volume fraction.