Estimating the discretization dependent accuracy of perfusion in coupled capillary flow measurements

One-compartment models are widely used to quantify hemodynamic parameters such as perfusion, blood volume and mean transit time. These parameters are routinely used for clinical diagnosis and monitoring of disease development and are thus of high relevance. However, it is known that common estimation techniques are discretization dependent and values can be erroneous. In this paper we present a new model that enables systematic quantification of discretization errors. Specifically, we introduce a continuous flow model for tracer propagation within the capillary tissue, used to evaluate state-of-the-art one-compartment models. We demonstrate that one-compartment models are capable of recovering perfusion accurately when applied to only one compartment, i.e. the whole region of interest. However, substantial overestimation of perfusion occurs when applied to fractions of a compartment. We further provide values of the estimated overestimation for various discretization levels, and also show that overestimation can be observed in real-life applications. Common practice of using compartment models for fractions of tissue violates model assumptions and careful interpretation is needed when using the computed values for diagnosis and treatment planning.


Introduction
Quantitative measurements of hemodynamic medical parameters based on tracer kinetic modeling are widespread both in research and in clinical practice [1][2][3]. Perfusion maps, as well as other parameter maps arising from tracer kinetic modeling can be combined with anatomical information and have proven to be of particular value in e.g. stroke studies or localization of trauma. Among the physiological parameters obtainable from tracer kinetic modeling, perfusion has been found particularly difficult to describe reliably on a voxel-basis [4]. These limitations are caused by issues in the numerical implementation [4], but might also depend a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 ordinary differential equation C 0 ðtÞ ¼ P a c a ðtÞ À P v c v ðtÞ; Cð0Þ ¼ 0: Here, c a , c v are the plasma CA concentrations at the inlet and outlet of O i and P a , P v is the normalized flow [mm 3 s −1 mm −3 ] at these locations. In the following, it is assumed that the plasma tracer concentration c a at the inlet is known. In clinical practice this can be accounted for by measuring c a in a feeding artery [9]. Since c v is usually unknown, additional assumptions need to be made if one wants to reconstruct the perfusion P a from a given tissue curve C. The convolution model and the maximum slope (MS) model diverge in further assumptions.

The convolution model.
For derivation of the deconvolution model, one approach is to assume there is an unknown probability distribution of transit times h through O i , cf. [1]. This leads to Combining this with (1) yields C 0 (t) = P a c a (t) − P a (h Ã c a )(t). Integrating and using basic properties of the convolution one obtains the general solution where the impulse response function I is defined as IðtÞ :¼ P a ð1 À R t 0 hðsÞdsÞ. The task of identifying I(t) given a tissue curve C(t) and an arterial input function c a (t) is a deconvolution problem. If I(t) is recovered, P a can subsequently be estimated as P a = max t I(t). There are several methods to perform the deconvolution. A standard approach using Fourier-based algorithms is sensitive to the presence of noise [9]. Another class of deconvolution algorithms gaining increasing attention are based on Bayesian modeling [10]. However the numerical handling is still difficult since complex and error-prone numerical integration has to be performed [10]. A popular class among deconvolution algorithms is based on singular value decomposition (SVD) [9]. These algorithms have shown to be robust for a reasonable noise level. Also, they can be easily adapted to be robust against delays in tracer arrival using block-circular structures (bSVD cf. [11]). In order to identify the impulse response function I(t) from applied data, we hence decided to use the bSVD model as proposed in [11].
2.1.2 The maximum slope model. In the MS model it is assumed that when c a has its maximum, only a negligible amount of CA is leaving the system [12]. For this time interval, (1) reduces to C 0 ðtÞ ¼ P a c a ðtÞ; Cð0Þ ¼ 0: One can see that if c a has a maximum, also C 0 must have a maximum since stationarity in P a is assumed. Hence, it holds that 2.2 A synthetic model for capillary flow model system we instead expect a set of coupled equations where each voxel can be regarded as an inlet for surrounding voxels. Hence, in order to make a realistic synthetic model for capillary flow, we decided to describe the CA propagation as a spatially coupled transport process, i.e. using partial differential equations (PDE) for transport. This PDE model is used for validation of the traditional models.
A major difference between our coupled flow model and traditional tracer kinetic modeling is the normalization of the flow field. To avoid a discretization dependent flow field for the PDE model, we instead of perfusion use vector valued surface fluid flux q = q(x) [mm 3 s −1 mm −2 ] as the fluid carrying quantity, in agreement with geoscience and porous media simulation theory. The fluid flux is a vector field describing the volume of fluid per unit time flowing across a sliced unit area of the sample. A detailed outline of how the flow field was obtained can be found in Section 2.7.

A model for indicator dilution
This section describes a model for CA propagation in the tissue. We assume that homogeneously dissolved CA is entering the domain along with the fluid flowing into the ROI via the source, and similarly extracted at a sink. In order to define meaningful and continuous contrast agent concentrations, we first describe CA concentration in an (arbitrarily) small tissue volume O ε . Let V ε be the volume of O ε centered around x and let v ε be the blood volume within the same control region. Letting the control region go towards zero volume, the porosity ϕ(x) ≔ lim ε!0 v ε /V ε [mm 3 mm −3 ] reflects the local, relative volume of the vascular system. The simplification of porosity as a continuous function is frequently performed in flow simulations. The flux q(x) as well as the porosity ϕ(x) are assumed to be stationary and hence independent of time. We further introduce C = C(x, t) and c = c(x, t) as the average CA concentration within V ε and v ε , respectively. By definition, we obtain the relation C(x, t) = ϕ (x)c(x, t). The rate of change of tracer molecules within a control volume O i can hence be phrased as d dt where the assumption of stationary ϕ(x) was used. Assuming mainly transport and marginal diffusion, the change in tracer mass within O i occurs from advective flow and the source and sink field Q(x). Let us write the source-and the sink term as Q(x) = Q si (x) + Q so (x) where Q si (x) < 0 is the sink and Q so (x) > 0 is the source, and zero elsewhere. Note that R O Qdx = 0. Using (6) and following the principle of conservation of tracer molecules, the rate of change of contrast agent in a control volume O i is modelled as Z where n(x) is the outward unit normal on @O i . Eq (7) is consistent with the continuity equation on local form for the initial condition c(x, 0) = 0 in accordance with no initial tracer at t = 0. Eq (8) is a linear transport equation in c(x, t). Following [13], (8) admits a unique local solution.

Relating the transport equation model with the traditional deconvolution model for perfusion
In this section we describe how the continuous model is related to the traditional deconvolution model. We will show that in the continuous model the flow into each voxel can be described as a traditional model with arterial input determined by adjacent upstream voxels. Let us start by modeling the CA concentration in a given voxel O i using traditional models. For sake of simplicity we assume that Q so = Q si = 0 within that voxel. It is possible to extend the following approach also to voxels where this is not the case. Define the outward normal vector n and voxel face areas of inflow and outflow over the boundary as S in ≔ {x 2 @O i : q Á n < 0} and S out ≔ {x 2 @O i : q Á n > 0} respectively. For the domain O i we define the arterial input c in as the weighted average of the tracer flux across S in c in ðtÞ :¼ We define local perfusion P v within O i as the total feeding fluid inflow divided by the volume, Given incompressible flow, the rate of fluid entering the region is the same as the rate of fluid leaving it. Further, let c i (t) denote the average fluid CA concentration within O i . Then it holds that P v = P out and we can describe c i by the traditional model (1) with solution The arterial input c in is determined by (9), which recursively depends on all upstream voxels until the global arterial input is reached. To verify this relationship numerically, we simulated a tissue curve C i (t) using both the continuous PDE model as well as analytical recursive convolution by (12). We refer to the latter approach as local convolution. The two curves have an almost perfect match, as seen in Fig 1 (left). It follows by recursion that the concentration at voxel i can be written as a convolution of the (global) arterial input function with a weighted average of all upstream impulse response functions. Deconvolving a tissue concentration C i with the global AIF will yield an impulse response function which depends not only on the local flow and porosity, but on flow and porosity of all upstream voxels. This relationship was also confirmed experimentally: Fig 1  (right) shows the impulse response function determined by analytical recursive convolution and the numerically achieved impulse response function obtained from deconvolving a tissue curve of the continuous model with the global arterial input function. The simulation was performed at location (1,20) of the digital phantom. The two curves coincide almost perfectly, highlighting the validity of the established theory.
These results show that the PDE model and the convolution model are equal in terms of local, voxelwise flow estimates if the convolution model is applied with the local arterial input. Also, the impulse response function obtained by convolution of the global arterial input function is identical to an analytical recursive convolution along all upstream voxels. This demonstrates that the perfusion recovered by traditional models depends on all upstream flow. However, for meaningful interpretation of the perfusion the entire streamline length within the capillary system needs to be taken into consideration. These results demonstrate that perfusion is depending on the geometry of the streamlines, hence the geometry of the capillary system.

Relating flux with perfusion
The flow model described in (21) uniquely determines the flux field q(x). However, in pharmacokinetic modeling the parameter of interest is usually CBF, which we will denote by P(x) as the voxelwise field of perfusion. Surface flux and perfusion are physically distinct, and there are at least two differences between q(x) and P(x). First, flux is a vector field and perfusion is a scalar field. Second, the flux is normalized to a surface area and the perfusion is normalized to a volume. Hence the flux describes flow over a surface separating spatial regions, while the perfusion describes blood leaving/entering a compartment within a given volume. According to the common understanding of perfusion, P(x) is the amount of blood feeding a tissue volume per unit time, with units [mm 3 s −1 mm −3 ]. In this work we address the fine scale setting, where the perfusion is taking place on a voxel level. At this level, a clearer understanding of how perfusion relates to the flux is desirable.
One straightforward approach for converting flux into perfusion could be to estimate the perfusion as the total inflow (or outflow) of fluid (e.g. arterial blood) into a control region per unit time, and then normalizing with the control region volume. This is a valid approach only if the control regions are not feeding each other, and is therefore well-founded for the entire organ, in line with the theoretical foundation of traditional compartment models for perfusion where a control region has its own source of feeding arterial blood, independent of neighbour regions.
On the other hand, if the control region is a single voxel or a sub-division of a capillary system with sequentially feeding arterial blood, the traditional model assumptions are violated since every control region will feed its neighbours, thus becoming a coupled system of flow. Simply summing the total inflow into a voxel and dividing by the voxel volume will overestimate the perfusion as the normalization refers to the wrong volume. This phenomenon is  (12) with experimental value of P v = 5328 ml/min/100ml at the given location and c in taken locally from upstream voxels around the simulated voxel. The two curves have an almost perfect overlap. Note that the numbers used for the perfusion is unrealisticely high since normalization is performed with respect to the volume of only one voxel. Right: Red curve shows the computed impulse response functions (IR) at location [1,20] using the global arterial input function. Blue curve shows the analytic impulse response function given by a convolution over all upstream flow. The two curves have an almost perfect overlap. These numerical experiments support that the computed impulse response function by traditional methods is not the directly feeding impulse response function, but rather a recursive impulse response function depending on all upstream voxels.  3 ]. However, for another discretization shown in the middle, the perfusion within each of these sub-volumes becomes P 2 = F 0 /V = 2P 1 . Taking the average across the two sub-volumes, it is clear that the perfusion is overestimated with a factor of two. A discretization dependent perfusion estimate is not recommendable, and the perfusion estimate of P 2 is clearly wrong.
In the following we introduce a meaningful notion of perfusion for the fine scale continuous model. To do this, we will consider distribution volumes which are following the streamlines, as shown in Fig 2 (right). For each point of a streamline we will select a small perpendicular disk with radius chosen in such a way that the total flow over each disk is constant along the streamline.
More precisely, let us consider an arbitrary streamline S O R 3 of length l > 0 and parametrization s: [0, l] ! S. We start by calculating the total flow over a small 2D disc perpendicular to the streamline. Let y 2 S be an arbitrary location along the streamline. The total flow F over a 2D disc B(y, R(y)) perpendicular to the flow field q(y), centered at y and with radius R : S ! R þ , is given by In order to calculate the perfusion, we need to establish the volume of a small tube around the streamline. We will not consider a tube with constant radius, but one with spatially varying Perfusion within a small volume. Left: A compartment with volume 2V is exposed to a flow F 0 [mm 3 s −1 ] of fluid. By definition, the perfusion within this compartment becomes P 1 = F 0 /(2V). Middle: The same volume is divided into two compartments (e.g. voxels), and the perfusion for each of the compartments becomes P 2 = F 0 /V = 2P 1 . Discrepancy between the two discretizations occurs because the flow is counted twice as it is fed from one voxel to the other. Right: As a solution to the described problem we rather pick out a true distribution volume ΔV (area in this 2D sketch), which is a small area around a given streamline along the centre line of the grey area. This is the true distribution volume (area in this 2D sketch) which is fed with arterial blood from the incoming fractional flow ΔF 0 . The correct perfusion within ΔV is therefore ΔF 0 /ΔV. The entire compartment can further be divided into similar infinitesimal distribution volumes, thus providing locally correct perfusion estimates. radii r : ½0; l ! R þ . The total volume of such a tube is given by Note that R(y) ≔ r(u) for some u 2 [0, l]. We define the perfusion at the arbitrary point y on the streamline as P s ðyÞ :¼ lim Fðy; εRðyÞÞ VðεrÞ for RðyÞ :¼ 1= In this expression the radii R(y) are chosen in such a way that in the limit when ε ! 0, the perfusion is constant along the streamline. To see this, let us assume that q is differentiable with Jacobian J. Using a Taylor expansion of q(x) around y, the Lagrange remainder theorem, as well as a change of coordinates z = (x − y)/(εR) yields where z i 2 (0, z i ) for every vector element i, and simplifications are due to RðyÞ ¼ 1= ffiffiffiffiffiffiffiffiffiffi jqðyÞ p j and n ≔ q(y)/|q(y)|. Combining this result with (15) yields what we refer to as global perfusion Note that (17) is independent of the spatial location y along the streamline, and is an explicit formula for converting flux into perfusion, showing that the perfusion scales with the streamline length l, as well as with the geometry of the domain, represented by the radii r(u).

A method to estimate local porosity
Porosity and CBV have the same definition, and we can therefore state that ϕ CBV. It is known from literature on traditional models [1] for perfusion that CBV for the entire compartment can be expressed as It is not obvious that (18) is valid also for a 1C field model where the voxels are feeding each other. We will now show that this is indeed the case. Let us switch to a discrete setting. Consistent with the considerations in Section 2.4, the CA concentration in any voxel can be described by C i (t) = ϕ i (J i Ã c in,i )(t), where the local arterial input is given by c in,i (t) = 1/P(P 0 c a (t) + ∑ j2J P j c j (t)) and J i (t) = (P/ϕ)e −(P/ϕ)t . In c in,i (t), J is the index set of all adjacent, upstream voxels and P = P 0 + ∑ j2J P j . Here P j is the normalized volume flow across voxel-face j [mm 3 s −1 mm −3 ] and P 0 > 0 if voxel i has arterial contribution. Furthermore, let us assume that q is a uni-directional flow field across each voxel face.
We will now use induction to show that R 1 0 c i ðsÞds ¼ R 1 0 c a ðsÞds, and then (18) follows. Let I k denote the set of voxels which have k layers of upstream voxels. E.g. I 0 is the set of all voxels, which have no upstream voxels, I 1 is the set of voxels which are fed by I 0 and so on. As an assumption, the same voxel can not be member of several I k , thus there is no flow interaction between voxels within the same I k . Induction will be carried out over k.
Induction basis: Let k = 0 and let i 2 I 0 be arbitrary. Since the area under the convolution of two functions equals the product of the area of its factors, R 1 0 c i ðsÞds ¼ R 1 0 c a ðsÞds and the claim follows.
Induction step: Consistent with our assumptions, for any voxel at location i 2 I k+1 which has the voxels J I k as their upstream neighbors, we find the following expression: Splitting the convolution integrals into separate factors, applying R 1 0 J i ðsÞds ¼ 1, as well as the definition of P yields the claim.

Simulation of capillary flow
In this section we describe how we simulate a flux field q(x) driving the transport of fluid and tracer. The modeling is in agreement with previous work on capillary perfusion simulations [14][15][16].
For the time being we will not consider contrast agent concentrations, but only the fluid flow in general. In-line with standard theory for a steady-state flow of an incompressible fluid and with Darcy's law [17], we assume that the fluid-flow q(x) obeys the following set of local PDEs Here Q [mm 3 s −1 mm −3 ] is the user-defined source-and sink term, which we assume to be only non-zero within the source or the sink, k = k(x) is the intrinsic permeability tensor, p(x) is the pressure, and μ(x) is the viscosity of the fluid. For simplicity, we will assume that k is a symmetric and positive definite tensor with only nonzero diagonal elements k ii = k, in accordance with a homogeneous porous medium. Using (20) yields the following elliptic partial differential equation in the pressure field p within the closed domain O, for the outward unit normal n(x). After solving (21) for the pressure p, the flux field was estimated according to (20).

Numerical implementation
First solving (21), and then (8), we set up a forward simulation of blood flow and indicator dilution through the capillary system. A standard arterial input function was chosen [9], the gamma-variate function c a (t) ≔ D 0 (t − t 0 ) α e −(t − t 0 )/β for default parameters α = 3, D 0 = 1 mmol l −1 s −1 , β = 1.5 s and t 0 = 0 s. Average, ground truth perfusion was chosen as 50 ml/min/100ml, a typical value for human brain perfusion [18,19]. The field of view was chosen as 2mm × 2mm × 2mm, in the order of the capillary bed or individual capillaries, ranging from 0.1 mm to 3 mm [14], or 0.25 mm to 0.850 mm [20]. The source term was assigned to the upper left voxel and the sink term was assigned to the lower right voxel. The source can be understood as the arterial compartment, the sink as the venous compartment, and the remaining field of view as the capillary system. The arterial input function (AIF) was measured in the source. Permeability was chosen to be isotropic and constant throughout the domain k = kI for the identity I and k = 5 × 10 −6 mm 2 . Dynamic blood viscosity was chosen as μ = 5 × 10 −6 kPa s according to [21]. Porosity (e.g. CBV) was assumed to be ϕ = 0.05, in line with measured CBV of the brain [19]. Eq (21) was solved numerically using two-point flux-approximation (TPFA), well known within porous media simulations [22]. The transport of CA described in (7) was implemented using first order upwinding [23], yielding a discrete 2D+time CA concentration map C(x i , t j ). From the porous media model using (21) and (8), streamlines to compute global perfusion P s were found from tracking of the flux vector field q by FACT [24], known from tractography.
Prior to reconstruction of perfusion using traditional models, the CA concentration map C (x i , t j ) was downsampled to a time-resolution of 0.1 s. In order to simulate different spatial resolutions of the scanning process, the data was averaged into blocks of {1, 2, 4, 8, 16, 32, 64} voxels with corresponding voxel sizes. Success of restoration was measured in terms of absolute, relative reconstruction error, RE(a, b) ≔ 100% Á |a − b|/b, where a is reconstructed values and b is ground truth values. Local perfusion P v was computed according to (10). Global perfusion P s was computed according to the streamline definition (17).

Reconstruction of perfusion within real data
In order to illustrate the effect of overestimation also in real data we applied the deconvolution model to a clinically acquired human perfusion CT dataset of a 56 years old male admitted with suspicion of stroke to the Radboud University Medical Center in Nijmegen, the Netherlands. The perfusion scan was obtained using a Toshiba Aquilon ONE scanner, voxel-size 0.43 mm × 0.43 mm, slice thickness 0.5 mm, contrast agent 50 ml Xentix 300, total scan-time 114 s, time resolution ranging from 2.1 s in the early-to 30 s in the late phase of CA uptake. Motion correction was performed with respect to the first timepoint using Euler transformations [25]. The arterial input function was manually selected within the middle cerebral artery (MCA) by a medical expert. Since we expected to see local overestimation effects mainly for small voxel sizes, the data was processed at full resolution (512 × 512 × 320 voxels). To cope with noise, we applied gaussian smoothing with standard deviation of 1 voxel and window size [5,5,5]. Relative concentrations were estimated from the CT signal assuming a spatially independent proportionality constant. The brain tissue was segmented automatically and used as ROI for the perfusion analysis.

Reconstruction of perfusion within synthetic data
Tracer dilution in the flux-field was simulated and from the resulting intensity time curves we tested the convolution based traditional model (bSVD) (3) as well as the maximum-slope (MS) model (5) for their capability to recover perfusion. Recovered perfusion maps P bSVD and P MS were compared against the two ground truth perfusion maps P s and P v depicted in Fig 3. As an internal control of P s , the average P s at maximal resolution was found to be 49.59 ml/min/ 100ml, for all practical means identical to the global input perfusion of 50 ml/min/100ml mediated through the source. Results from reconstruction of the porosity ϕ (i.e. CBV) according to (18) resulted in reconstruction errors of <1% for all voxel sizes.
We performed two different normalizations of the restored flow, a normalization (i) with respect to volume, and (ii) with respect to surface. The volume normalization (i) implies normalizing the flow to [ml/min/100ml], most in line with common units for perfusion. A comparison of ground truth perfusion to reconstructed perfusion using volume normalization is shown in Fig 4. For a voxel size corresponding to the entire ROI (voxel size = 3 mm) the reconstructed perfusion of P bSVD and P MS is close to the ground truth perfusion P s and P v . For any voxel size smaller than the entire domain the relative error increases inversely with voxel size, in particular for reconstruction by bSVD. Global perfusion P s is not depending on voxel size.
For surface normalization (ii) we first computed the absolute flow F [ml s −1 ] of ground truth perfusion as well as reconstructed perfusion, and then normalized the flow to the surface area of the distribution volume, F/S, here referred to as surface normalized flow. This interpretation of perfusion has no standard clinical unit, thus we choose to scale the surface area S to [mm 2 ]. Reconstruction results of surface normalized flow is shown in Fig 5. For a voxel size corresponding to the entire ROI (voxel size = 3 mm isotropic) the surface normalized flow of P bSVD and P MS is close to the ground truth. For any voxel size smaller than the entire domain the relative error increases inversely with voxel size, in particular for reconstruction by bSVD. Both global perfusion P s and local perfusion P v are dependent on voxel size in the framework of surface normalized flow.

Reconstruction of perfusion within real data
Perfusion for the entire brain by averaging the concentration time curves first and then performing the bSVD yielded a perfusion of P bSVD = 24.79 ml/min/100ml. As a second step, voxelwise perfusion was estimated, depicted in Fig 6. These values yielded an average perfusion of " P bSVD ¼ 64.36 ml/min/100ml, corresponding to an overestimation of perfusion with RE = 159.60% compared to the value obtained for the entire brain.

Discussion
It has previously been shown that perfusion reconstructed from traditional 1C models in a coupled system is discretization dependent (cfr. Fig 2) [6][7][8]. As a consequence, the obtained results will strongly depend on acquisitions parameters and post-processing tools. It is unknown to which extent the pharmacokinetic modelling overestimates perfusion and whether the error is homogeneously distributed or not. Considering this, the shortcoming of existing perfusion formulations has not been sufficiently well accounted for within clinical studies [26,27]. To clarify the potential impact of limitations seen within existing perfusion models, our main contribution in the current work is to quantify the observed error. To the best of our knowledge, such quantification has not been carried out previously.
Our results strongly support the usage of traditional 1C models for entire regions exclusively fed by the measured arterial input. Moreover, our results also show that when traditional models are applied only to parts of the system, the measured perfusion is overestimated (cfr.  Observed error in perfusion for a voxel size of *2 mm was found to be * 40% for reconstruction by bSVD and * 20% for reconstruction by MS (cfr. Fig 4). This is a relevant spatial scale for today's MR acquisitions. The error is expected to increase in future acquisitions along with hardware and software improvements leading to higher spatial sampling.
There are at least two reasons for overestimation of perfusion in traditional 1C models. The first reason is that blood passing through a voxel without being locally delivered to the capillary tissue will contribute to artificially high perfusion values. This issue has not been accounted for in our digital model but my be handled by more complex multi compartment spatial modelling described in i.e. [8]. The second reason is thoroughly described here, and relates to estimation of an incorrect distribution volume used for computing the perfusion. Overestimation of perfusion obtained within the digital phantom was also confirmed by real data experiments, where we showed local overestimation of perfusion for voxelwise estimates as compared to an averaging of concentrations for the entire volume of interest (cfr. Section 3.2).
In order to demonstrate our results we introduced two definitions of voxelwise perfusion, global perfusion P s and local perfusion P v . Local perfusion P v is in line with [7] where the authors demonstrated a discretization dependent flow without connecting it mathematically to perfusion. Theory and examples in our work show that this definition of perfusion does not comply with the physical understanding of perfusion as a feeding arterial blood flow. The correct distribution volume is not accounted for and the obtained perfusion will be strongly overestimated compared to the actual perfusion. However, our analyses show that traditional models would restore the local flow value if the local arterial input function was selected, implying that traditional models are accurate as long as the model assumptions are not violated. The coupling between the continuous porous media model and the convolution model in Section 2.4 demonstrates that there is no contradiction between these two models. The problematic issue of traditional models is related to physical interpretation and normalization with respect to incorrect distribution volume.
Global perfusion P s models perfusion along the streamlines and most accurately reflects the physical notion of volume flow within the correct distribution volume according to mathematical definitions. We showed that P s is independent of discretization (cfr. Fig 4), P s is a constant quantity along the streamline, and scales with streamline length and geometry according to (17). For our purpose, the concept of P s was useful as a realistic ground truth in order to clarify the definition of perfusion as a flow that must be normalized along the entire capillary length, where the blood undergoes a transition from arterial to venous blood. Traditional perfusion measurements are thereby dependent on both the discretization and the geometry of the domain. While the geometry in our experiments is limited to a simplified square with diagonal flow gives fairly simple analysis, estimation of P s in real applications is practically difficult due to a highly complex and unknown microvasculature. Nevertheless, the derivation of the theoretical relation between P s and streamline length is valid for an arbitrary geometry.
Development of new field models for perfusion is highly demanded due to the scale dependency, and a few initiatives have been proposed based on multi compartment flow models [8,16]. Alternative approaches with estimation of voxel-speciffic AIFs have also been suggested [28,29]. These approaches may suffer from instabilities or errors which have not yet been handled [29]. A possible hybrid approach that may have clinical relevance is to estimate a modest number of local AIFs feeding anatomically sensible regions on a scale way above the capillary systems. This will however require some prior knowledge of local anatomy and vasculature.
It was previously suggested to normalize the flow by surface instead of volume [7]. Our experiments suggest that a surface normalization is nevertheless discretization dependent, and traditional 1C models are not able to restore this type of perfusion, neither for global, nor for local perfusion (cfr. Fig 5, blue and black lines).
We have also shown there is a global error directly scaling with smaller voxel sizes (cfr. Fig 4). A comparison between individual scans with otherwise equal acquisition parameters and post-processing chains should ideally adjust for the global error in the interpretation of absolute perfusion values. As such, voxelwise maps of perfusion could still be of high clinical value as the main goal is a comparison of perfusion between patients or between repeated scans of individuals. However, particular care should be undertaken in the case of comparing perfusion data of various resolution. Multicenter or retrospective studies are particularly susceptible to this issue where data collected from various sources are different with respect to hardware, resolution and post-processing tools that can affect the discretization level. Future study design should account for this limitation and special care should be undertaken to ensure equal level of discretization in perfusion estimates.
In addition to the observed global error we also observed inhomogeneous reconstruction errors within the capillary system. This becomes clear in Fig 3 where the reconstructed perfusion maps P bSVD and P MS are strongly unlike the ground truth perfusion maps of P s and P v . This inhomogeneity leads to locally inaccurate estimates of perfusion within a patch, even if the global average value within the patch were correct. If analyses of voxelwise perfusion is undertaken with high resolution the local error can become large within each capillary patch.
Regarding the CBV estimates, estimation of blood volume is stable, and varying voxel size had little impact on the results. These results are in agreement with the analyses in Section 2.6 supporting the usage of (18) for computing voxelwise CBV with high accuracy for any voxel size.

Conclusion
Our experiments confirm that traditional 1C models for perfusion perform well if they are applied to the entire domain. However, when they are applied to fractions within a coupled domain, perfusion becomes scale dependent. We quantified substantial and increasing reconstruction errors of perfusion as a function of smaller voxel size, and we also found similar effects in real data. The observed reconstruction error for a simplified geometry in clinically relevant resolution was between *20% to *40%. The error is expected to vary with geometrical complexity of the capillary system and increase with higher spatial resolution. The reason for the observed errors is not numerical instabilities in the deconvolution but rather that traditional 1C models will not account for correct distribution volume within smaller fractions of the region of interest. As a consequence, interpretation of absolute perfusion for the purpose of diagnosis and disease monitoring must be undertaken with care, and comparison of perfusion values from individual dynamic data sets with different resolution is not recommended.