A calibrated physical flow standard for medical perfusion imaging

In the medical sector, various imaging methodologies or modalities (e.g. MRI, PET, CT) are used to assess the health of various parts of the bodies of patients. One such investigation is the blood flow or perfusion of the heart muscle, expressed as the (blood) flow rate normalized by the mass of the volume of interest. Currently there is no physical flow standard for the validation of quantitative perfusion measurements. This need has been addressed in the EMPIR 15HLT05 PerfusImaging project. A phantom simulating the heart muscle has been developed with the capability that it can reproducibly generate a flow profile with individual flow rates known with a relative uncertainty of about 10% ( k = 2) and total flow rate known with an uncertainty of 1% ( k = 2). An overview of the phantom and its validation is given. Next, a new analysis method is presented to analyse the sequence of images which are acquired when using a standard dynamic imaging protocol. It is concluded that the new, alternative approach gives results comparable to the standard analysis method.


Introduction
In the medical sector, various imaging methodologies or modalities (e.g. MRI, PET, CT) are used to assess the health of various parts of the bodies of patients. One such investigation is the blood flow or perfusion of the heart muscle, expressed as the (blood) flow rate normalized (divided) by the mass of the volume of interest, with unit mL/min/g. A decreased perfusion of the heart muscle or myocardium is an imaging biomarker for increased risk for a heart attack and can cause chest pain. The current state of the art is that MRI images are qualitatively judged by physicians. Quantitative perfusion values are not used and uncertainty evaluation of these values does not yet exist.
There is no physical flow standard for the assessment and validation of myocardial perfusion imaging methodologies, resulting in a large proportion of medical diagnoses being inaccurate and highly dependent on the scanner type, software used and the clinical operator. In the EMPIR 15HLT05 PerfusImaging project [1] a phantom simulating myocardial perfusion has recently been developed with which imaging modalities can be tested. The design was based on an earlier, already existing phantom made by KCL [2].
In this paper the construction and validation of the phantom is described which involved several iterations with design updates, computational fluid dynamics (CFD) simulations, 3D printing of the phantom, ultrasound imaging velocimetry (UIV) and MRI measurements involving all partners. Dynamic contrast enhanced MRI (DCE-MRI) measurements were performed as well. The standard and a new flow model and analysis method to interpret the imaging data will be presented.

Phantom design
The phantom [3] mimics the complete flow through a human heart including the four heart chambers, pulmonary volume and the blood flow through the heart muscle. This design has been chosen in order to mimic the clinical practice as faithfully as possible (note that the contrast agent is injected upstream of the heart). The part simulating the tissue of the heart muscle or the myocardium is of special interest in this study. The arteries in this tissue are simulated by more than 200 parallel channels or capillaries with cross-sectional area ranging from 1 to 7 mm 2 . This configuration mimics the multitude of arteries in human tissue which have different sizes and different flow resistances. Upstream of the channels is a pre-chamber with a volume of 36 mL which is meant to enhance equal repartition of the flow and contrast agent to all channels. This myocardium part was produced by 3D printing and is shown in Fig. 1. The engineering and production of the phantom was done by the company Zurich MedTech (ZMT). In Fig. 2 a schematic drawing of the phantom is shown including the location of the MRI imaging plane. The flow rates through the aorta (cardiac output) can be varied between 2 and 5 L/min and the perfusion rate through the myocardium between 1 and 5 mL/min/g, corresponding to flow rates between 85.3 and 426.5 mL/min. In this study, the flow rate through the aorta was set around 4 L/min and the full range of flow rates through the two myocardia were tested.
The design combines the aspect of being a myocardium phantom mimicking the human heart with the aspect of being a (somewhat special type of) flow standard by generating a reproducible flow pattern with known flow rates through the myocardium channels. The total flow rate follows from the reading of a calibrated flow meter, and the repartition amongst the channels follows from PC-MRI measurements described in section 2.3 and supported by CFD simulations described in section 2.1. In the next sections we will rather focus on relative flow rates, because results are then in theory independent of the absolute flow rate and relative values have a greater clinical interest than absolute values. The difference to already existing phantoms is exactly this variation in flow rates over the cross section induced by length and cross-sectional area differences between the capillaries.

CFD simulations
The design goal was a channel configuration with a rotationally symmetric flow profile with maximum flow velocity at the centre and linearly decreasing flow velocity to − 30% at the outer wall (for channels with identical cross sections). While Poiseuille's law for a laminar flow regime suggests a linear increase of +30% in channel length towards the outer wall to obtain the desired linear variation in flow velocity, the OpenFOAM CFD solver was used by VSL to confirm this expectation and to ensure that the flow regime would be rotationally symmetric, at least in the model. The final design had four inlets from the sides to the prechamber in which a laminar flow element upstream of the channels was installed to remove any form of potential pressure gradient and turbulence in the cross section (the Reynolds number is below 250 in the prechamber). In Fig. 3 a plots of the results of the CFD simulation is shown. Approximately 4 million cells were used. For the large octagonal channels the CFD results are in line with the simple analytical Poiseuille model, see Table 1. For the smaller square channels with side lengths 1.0, 1.2, 1.4 and 1.8 mm, the predicted flow rates by CFD are smaller than the analytical ones, which is probably due to the mesh size of the CFD simulation.

UIV measurements
In order to measure the actual flow rates through the channels of the phantom, measurements with ultrasound imaging velocimetry (UIV) were performed by TU Delft. The earlier version of the phantom was based on thin walled straws [2]. The current version was 3D printed and had thicker walls [3], which turned out to be problematic for the UIV measurements. Some indicative results were finally obtained by inverting the direction of the channel bundle inside the compartment. The results suggest linearly decreasing flow velocities with radius. However, the imaging quality and reproducibility of the UIV measurements were so poor, that the target measurement uncertainty of 10% (k = 2) could certainly not be met using this technique and a detailed uncertainty analysis for the UIV measurements has therefore not been performed. In Fig. 4 a schematic view of the measurement set-up and a    Table 1. measurement image are shown.

Phase contrast MRI measurements
As alternative to the UIV measurements phase contrast (PC) MRI measurements were performed with an estimated expanded relative uncertainty of 10% (k = 2) for the flow velocities. This uncertainty has been derived from at the one hand comparing the values of the integrated flow velocity fields (using the known pixel size) with the readings of the calibrated flow meter, which has an expanded uncertainty of 1% (k = 2). At the other hand studies like [4] suggest a typical uncertainty of 10% for the full velocity field as measured by PC-MRI. There has been no traceable calibration of the full velocity field of the particular MRI device used for reasons of required effort and budget. The fact that the measured velocity field of the phantom agrees with the results predicted by CFD within 10% further supports the uncertainty claim of 10% for the PC-MRI device (and the correct working of the phantom).
Five flow rates between 85.3 and 426.5 mL/min were each measured twice, and the results were analysed in various ways by KCL and VSL. In Fig. 5 one such analysis by VSL is shown. The channel geometry is slightly different than that of Fig. 3, but the expected decrease in flow rate is the same. The octagonal channels on a horizontal line were identified by an algorithm, and the flow rates through channels were calculated by integrating the velocity field and expressed as a percentage relative to the flow rate in channel 5 so that the results for different flow rates can be compared. A decrease of 30 ± 10% with respect to the maximum flow rate is visible for all flow rates, be it not completely symmetrical between left and right side.
In summary, the reasoning regarding the 10% (k = 2) uncertainty and 'traceability' of the flow rates per channel of the phantom flow standard is as follows (note that the standard is not traceable in a strict metrological sense): • The calibration of the flow rates in the (large) channels of the phantom is done by PC-MRI measurements. • PC-MRI velocity field measurements are generally believed to have an uncertainty of about 10% (k = 2) (see e.g. Ref. [4]), and this estimate is supported for the particular PC-MRI device we used by the fact that the total flow rate as measured by a calibrated flow meter (U = 1% (k = 2)) agreed within 10% with the integrated velocity field. • Analytical results and results predicted by CFD are in line with what has been measured within 10%, supporting our confidence in PC-MRI and the phantom.   • The spread of 10% visible in Fig. 5 is mainly due to comparing the normalized results for different flow rates. The spread for individual (absolute) flow rates is generally smaller than 3.2% and therefore doesn't increase the rounded uncertainty value of 10%.
This expanded uncertainty is largely fit for purpose when comparing the phantom reference values with values obtained using a clinical test protocol involving DCE-MRI (see section 3) and therefore a more precise analysis was not judged necessary.

Dynamic contrast enhanced MRI protocol
For clinical measurements of myocardial perfusion in patients a socalled Dynamic Contrast Enhanced (DCE) MRI protocol is being used. A bolus of contrast agent (CA, e.g. gadolinium) is injected in the patient upstream of the heart and this is simulated by injecting such a bolus at a specific point in the phantom representing the vena cava. The CA concentration in the blood is roughly proportional to the measured MRI signal. The bolus of CA mixes with the blood (or water) and passes through the four heart chambers. Then, the MRI signal is measured as function of time in (a pipe simulating) the aorta, upstream of the myocardium. The resulting time series is called the arterial input function (AIF). In the same image, the concentration inside the myocardium is also visible. After conversion from signal intensity to CA concentration, the perfusion of the myocardium can be estimated. In patients typically three slices of the heart are being measured, though one slice at a time is used during the analysis phase. In the experimental work with the phantom one slice of the phantom has been measured and analysed. To reduce the effect of the non-linearity of the relationship between CA concentration and MRI signal intensity a dual bolus scheme is typically used. In this scheme a first 'pre-bolus' diluted in saline at a ratio 1:9 is injected from which the AIF is derived and scaled by a factor 10. Subsequently, a second undiluted bolus is injected from which the myocardium signal is taken. As the CA concentration in the myocardium is much lower (and more spread out over time) than in the aorta, the MRI signal magnitudes of AIF and myocardium measurements are similar, and the effects of non-linearity and signal saturation are diminished.
In the next sections the standard method and a new method for analysing the data are presented.

Standard method
The standard theory and method [5] for quantification of the perfusion assumes a system of volume V sys with one inlet and one outlet, see Fig. 6. The AIF is measured at the inlet by MRI-A (aorta) and is proportional to the CA concentration c in (t). The average CA concentration in V sys at time t equals c sys (t) and is measured by MRI-B, the measurement of the myocardium. The stationary flow rate through the system equals q in = q out and the quantity of interest is the perfusion f = q in /V sys .
Based on a mass balance of CA in the system, the following equation Assuming a linear and stationary system with impulse response function h(t), c in (t) = c sys (t) = 0 for t < 0 s and using q in = q out the outlet concentration can be written as where * denotes the convolution operation. Defining h(s)ds) and performing some transformations one can derive that Solving convolution Equation (3) for In practice a correction factor has to be applied for the accessible part of the myocardial tissue for the CA ('volume fraction'), or the 'nonplastic fraction' in the case of a phantom.
Advantages of the standard approach are its relative simplicity and the fact that the volume V sys doesn't need to be known. Possible limitations are that the measured concentration by MRI-B is in practice (at least for the phantom) rather an outflow CA concentration c out (t) than the system average concentration c sys (t), and interaction between different parts of the myocardium is not modelled.

Alternative method
An alternative model and method for calculating the perfusion is based on Fig. 7, where a system with a common pre-chamber is shown and two subsequent compartments in parallel. Note that this model is not limited to three compartments but can be generalized to an arbitrary number of compartments. For ease of presentation only three compartments are used in this paper.
It is again assumed that at location MRI-A the MRI signal is proportional to the CA inlet concentration c in (t). However, at MRI-B 1 and MRI-B 2 the signal is assumed to be proportional to the respective outlet concentrations. The perfusion f i in compartment i is defined by f i = q i /V i , the impulse response function for compartment i is denoted by h i (t) and the mean transit time T i of a unit of CA through a compartment i is given by  measured MRI signals are which is similar to Equation (3). In order to be able to calculate the perfusions an additional hypothesis is needed. The assumption made is that there exists a constant tissue delay factor d for which following relation holds true for any compartment i: From Equation (5) then follows that f i = d/T i . Using q in = q 1 + q 2 all perfusions as well as flow rates can be calculated when the measurement results c in (t) and c i (t) are available.
When a dual bolus scheme is being used or when there is a significant delay before the CA arrives at the myocardium additionally an additional time offset τ has to be used in Equation (5), leading to T i = τ + d V i /q i . Note that this is also the case for the standard method. In the assessment of the methods discarding τ and only using the main bolus gave better results and these are shown in Table 2.
Advantages of this alternative method are that it seems to be more realistic and thus might yield more accurate results, as it models the measurements at MRI-B i as proportional to the compartment outflow concentration rather than the compartment's average concentration, and it models a common large blood vessel (or phantom pre-chamber) upstream of smaller myocardium blood vessels (or phantom channels) in parallel. Disadvantages are the required knowledge of the (relative) volumes of the compartments (which is possible for a phantom but very difficult for tissue of patients) and the hypothesis of constant tissue delay factor (of which the exact value is required if absolute values of perfusion are of interest). Another disadvantage is that determination of the mean transit time requires integration over the complete relevant time domain of the signals c in (t) and c i (t), which is possible for a phantom, but problematic for patients where CA may re-enter the aorta and myocardium resulting in an undesired 'second-pass' perfusion MRI signal. Another assumption is that c 0 (t) is homogeneous at the outlet plane of the initial compartment. (A similar assumption is used in the standard method.)

Comparison of data analysis methods
The two analysis methods described in section 3.2 and 3.3 were applied to DCE-MRI data acquired for the phantom described in section 2. The target quantity was the ratio r of perfusion through the outer ring to the perfusion through the inner circle, i.e. r = f out /f in , see also Fig. 8. The two segments have equal area.
The channels close to the outer boundary were taken out in this analysis, as measurement results were unreliable in this area due to partial volume effects.
The results of the calculations for both methods are shown in Table 2. In this case both methods turn out to perform equally well and are both in line with the reference value based on PC-MRI measurements (which are in-line with CFD and analytical calculations).

Discussion
In section 2 a myocardium phantom was presented which has known flow rates (U = 10% (k = 2)) through the channels and which constitutes in that sense a 'medical flow standard' with which imaging modalities involving clinical protocols like DCE-MRI can be tested. The correspondence between the values as derived by the DCE-MRI measurements and the reference perfusion values are generally good (see Ref. [3] for more details).
Once a sequence of images has been acquired, a significant amount of post-processing is needed in order to derive perfusion values. Various methods exist and have been compared in literature [6]. In section 3 the standard method [5] and a new alternative analysis method were presented. The two analysis methods both rely on the validity of some additional assumptions. One of these assumptions is that the CA concentration at the inlets of all channels are equal at the same point in time and that the images are rotationally symmetric. Inspection of the raw MRI images directly shows that this is not the case, see Fig. 9. At the bottom and towards the right upper corner the MRI signal appears to be higher. As the contrast agent has a higher density than the water some stratification seems to happen especially in the case of the lowest flow rates. Updated designs with passive mixers may improve the situation. Another practical issue that was encountered was that removing all the Table 2 Results of the calculations using the standard method (r A ) and the alternative method (r B ). The reference value based on PC-MRI measurements is r A = 0.87 ± 0.05.   air bubbles from the phantom was not so straightforward. The comparison of the performance of the methods may be obfuscated by effects such as stratification effects as discussed in the last paragraph. In the analysis algorithm various parameters need to be selected, such as an image base line value, the pixels to be included in the analysis and similarly for the selection of time points. Each of these choices has an effect on the results which in some cases is non-negligible. These issues are also encountered in medical practice and the phantom can help to assess the importance of each parameter.
The goal of the phantom flow standard is to make MRI perfusion measurements more quantitative and more comparable across different MRI manufacturers, software settings and operators when measuring real patients. For understanding the flow physics through the phantom better, the alternative method can be useful. However, even if the alternative method would turn out to outperform other methods in a next version of the phantom, its requirements in terms of necessary (relative) volume information and measurement time range and other assumptions may limit its applicability in a medical setting. Application to sample patient data is a necessary next step to assess the usefulness of the alternative method in a medical context.

Conclusion
In this paper a calibrated physical flow standard for medical perfusion imaging was presented, i.e. a phantom mimicking a human heart with known flow rates through its channels which can be used to assess the performance of MRI scanners as well as of PET and CT scanners.
Reference values for the flow rates through the channels were measured with phase contrast MRI with an uncertainty of approximately 10% and were in line with CFD and analytical calculations. UIV measurement results were only indicative due to the thickness of the 3Dprinted material. Overall flow rates are being measured with flow meters.
A standard and an alternative, new model and associated analysis methods were presented that can be used to relate dynamic contrast enhanced MRI measurements with the reference flow rates. The methods performed approximately equally well for the phantom measurement data. Current research in the project consists of measuring the phantom with different scanners, extending the phantom to a multi-compartment phantom simulating more faithfully human tissue and applying the alternative data analysis method to real patient data in order to assess its usefulness in a medical setting.
For the future, it's envisaged that the optimized phantom flow standard (commercialized by ZMT) and related data processing methods will be used for scanner testing by medical centres, which will be a big step towards introducing quantitative perfusion measurements which are comparable between centres and imaging modalities.