## Journal of the American College of Cardiology

# Accurate noninvasive quantitation of blood flow, cross-sectional lumen vessel area and wall shear stress by three-dimensional paraboloid modeling of magnetic resonance imaging velocity data

## Author + information

- Received November 11, 1997.
- Revision received March 26, 1998.
- Accepted April 9, 1998.
- Published online July 1, 1998.

## Author Information

- Sten Oyre, MD∗,†,§,* (skejsao{at}aau.dk),
- Steffen Ringgaard, MS†,
- Sebastian Kozerke, MS§,
- William P. Paaske, MD, DMsc∗,
- Mogens Erlandsen, MS‡,
- Peter Boesiger, MS, PhD§ and
- Erik M. Pedersen, MD, PhD∗,†

- ↵*Address for correspondence: Dr. Sten Oyre, Department of Cardiothoracic and Vascular Surgery T, Aarhus University Hospital, Skejby Sygehus, DK-8200 Aarhus N, Denmark

## Abstract

Objectives. We present a new method in which a priori knowledge of the blood velocity fields within the boundary layer at the vessel wall, combined with acquisition of high resolution magnetic resonance imaging (MRI) blood velocity data, allow exact modeling at the subpixel level.

Background. Methods are lacking for accurate, noninvasive estimation of blood flow, dynamic cross-sectional lumen vessel area and wall shear stress.

Methods. Using standard acquisition of MRI blood flow velocity data, we fitted all data points (n = 69) within the boundary layer of the velocity profile to a three-dimensional paraboloid, which enabled calculation of absolute volume blood flow, circumferential vessel wall position, lumen vessel area and wall shear stress. The method was tested in a 8.00 ± 0.01-mm diameter glass tube model and applied in vivo to the common carotid artery of seven volunteers.

Results. In vitro the lumen area was assessed with a mean error of 0.6%. The 95% confidence interval included the specified tube dimensions. Common carotid mean blood flow was 7.42 ml/s, and mean (standard error) diastolic/systolic vessel area was 33.25 (0.72 [2.2%])/43.46 (0.65 [1.5%]) mm^{2}. Mean/peak wall shear stress was 0.95 (0.04 [4.2%])/2.56 (0.08 [3.1%]) N/m^{2}.

Conclusions. We describe a new noninvasive method for highly accurate estimation of blood flow, cross-sectional lumen vessel area and wall shear stress. In vitro results and statistical analysis demonstrate the feasibility of the method, and the first in vivo results are comparable to published data.

Blood flow velocity can accurately be determined noninvasively and in vivo by magnetic resonance imaging (MRI) phase contrast velocity mapping techniques (1–3). In each image pixel, velocity information can be acquired across an imaging plane, and all three orthogonal components of the velocity vector can be determined (4,5).

Translation of velocity data provides volume blood flow (1–3), and attempts have been made to derive methods for determination of wall shear stress (6,7), vascular compliance (8–10), and blood pressure (10). Correct determination of these variables is dependent on exact determination of the circumferential lumen vessel position throughout the heart cycle, but a reliable method with subpixel resolution has not been developed (2,6–10).

We previously presented a method for determination of wall shear stress based on a simple edge detection approach (6). The present report introduces a new concept for determination of cross-sectional vessel area and wall shear stress based on circumferential subpixel edge detection. A priori knowledge of the parabolic blood velocity distribution within the thin boundary layer at the vessel wall (11,12) is applied to velocity data from the whole circumference of the vessel. By this technique, a three-dimensional paraboloid (3DP) model can be fitted to the large number of velocity data obtained within this parabolic boundary layer. Position of the lumen vessel wall, wall shear stress and volume blood flow can then be determined at a subpixel level from the 3DP model.

The present study sought to 1) assess the validity of our 3DP model by in vitro studies and statistical analysis; and 2) apply the technique to the common carotid arteries of normal subjects.

## Methods

### Theory of 3DP model

The method is based on the following assumptions: 1) Blood velocity at the vessel wall is zero; 2) the blood velocity profile in the boundary layer of the common carotid artery is parabolic and has rotational symmetry; and 3) arteries are circular in shape when cut perpendicular to the longitudinal axis. These assumptions result in the following equation describing the 3DP blood flow velocity profile in the boundary layer close to the artery wall:

### Definition of vessel

Initially, a simple semiautomatic definition of the vessel wall was performed: A circle was placed through two points positioned manually on the “top” and “bottom” of the vessel by visual inspection of the magnitude image, as seen in Figure 2a.

### 3DP fitting

The mean in vivo boundary layer thickness was measured to be 1.6 mm in peak systole. The pixels for the 3DP fitting were selected in the middle of the boundary layer at a distance of 0.3 to 1.3 mm from the vessel wall, as illustrated in Figure 2, b and c. The pixels from this selection were subsequently fitted by multiple linear regression (least squares method) to the 3DP equation.

### Edge detection and blood flow calculations

The circumferential vessel wall position was computed throughout the heart cycle assuming zero velocity at the vessel wall, and the cross-sectional lumen vessel area was calculated. Volume blood flow was calculated by summation of all pixel velocities within the vessel. The peak central blood flow velocity was found as the average of the central nine pixels, that is, nine pixels covering a 1.5 × 1.5-mm^{2} square area in the center of the vessel (Fig. 2b).

### Wall shear stress calculations

*Wall shear stress* (WSS) is the term for the mechanical stresses on the vessel wall exerted by the flowing blood: *Wall shear rate* was calculated by differentiating the 3DP model at the wall (Fig. 2c). The peak systolic, end-diastolic and average WSS in the cardiac cycle were calculated. The peak systole and end-diastole heart phases were identified as the heart phases with the highest and lowest flow, respectively. To allow comparison with other studies, WSS was also estimated using the central peak velocities and applying the Poiseuille concept (13,14).

### Statistical analysis

The 3DP model statistics used were the standard deviation (SD) and the adjusted coefficient of determination (R^{2}). To provide the asymptotic standard error (SE) for the center (h, k), radius r and wall shear rate, we described the paraboloid in terms of the four variables: center (h, k), radius r and wall shear rate. These variables were determined using nonlinear regression (Levenberg-Marquardt method) and were used to calculate the asymptotic SE for WSS and radius (15). The radius SE was used to extract the diameter and lumen area SE.

### Applications of 3DP model

The model was applied in vitro and in vivo.

#### In vitro

A 1-m long glass pipe (Schott Geräte, Mainz, Germany) with an internal diameter of 8.00 ± 0.01 mm (factory specifications) was used in a steady manganese chloride–doped water flow loop. The pipe was surrounded by doped water. Flow was measured using the stopwatch and graduated cylinder method. Measurements were performed at 4.37 ml/s (Reynold’s number [Re] 695) and 6.45 ml/s (Re 1,027). For each flow rate, we made 10 measurements during which the glass tube was moved to different in-plane positions. The MRI variables were identical to those of the in vivo measurements (see next section), except for the maximal velocity encoding of 30 cm/s.

#### In vivo

Seven healthy, young volunteers were studied (mean age 26.4 years, range 20 to 36; 4 men), as approved by the institutional committee on human research, and individual written informed consent was obtained according to the Helsinki II declaration. Measurements were performed using a 1.5-T whole-body scanner (Philips Gyroscan-ACS-NT, Philips Medical Systems, Best, The Netherlands). To optimize the signal to noise ratio, a standard 8-cm diameter circular surface coil was positioned over the carotid artery inside a standard head coil. Perpendicular blood flow velocity measurements in the common carotid artery (CCA) were performed 2 cm below the carotid bifurcation after initial visualization of morphology by sagittal, transverse and coronal scout images (16). A standard electrocardiographic gradient-echo pulse sequence with bipolar velocity encoding gradients was used with the following variables: Thirty-two frames were recorded throughout the cardiac cycle with an interval of 25 ms, depending on the heart rate using retrospective triggering. The thickness of the slice was 7 mm, and in-plane pixel resolution was 0.5 × 0.5 mm^{2} (data acquisition matrix 128 × 128 points, 64-mm field of view). Two signal averages, 8-ms echo time acquiring a full echo and a maximal velocity sensitivity of ±90 cm/s were used.

## Results

### In vitro results

The flow, cross-sectional area, diameter and WSS for the two different flow rates are presented in Table 1. It can be seen that the lumen area was assessed with a mean error of 0.6% (diameter error 25 μm [0.3%]); 78 pixels were used in all 3DP fits. The in vitro mean adjusted R^{2} and SD for all fits were 0.978 m/s (range 0.966 to 0.984) and 0.0035 m/s (range 0.0026 to 0.0047), respectively. Nonlinear statistical analysis revealed for both the low flow (diameter ±95% confidence limits 7.961 ± 0.068) and high flow (8.010 ± 0.048) used that the confidence interval included the factory-specified lumen dimensions. WSS had a mean SE of 0.0015 N/m^{2} (range 0.0011 to 0.0020).

### In vivo results

The flow profiles in peak systole were blunt in all cases (Fig. 3). In late systole, a small area of blood flow separation was seen in six subjects. During early diastole, the dicrotic notch was evident as a second peak in the velocity patterns of all subjects. At this time, all profiles were still blunt and slightly skewed. During diastole, the velocity profiles became almost developed and parabolic.

Figure 4 shows the two-dimensional velocity profiles from peak systole and end-diastole from one subject as well as the result of the 3DP fit of the velocity in the boundary layer presented in a two-dimensional profile overlaid with the measured pixel velocities.

The in vivo results for peak systole, end-diastole and mean values of all heart phases are shown in Table 2. Table 3summarizes the in vivo statistics from the 3DP fits and nonlinear statistical analysis.

Figure 5 shows blood flow and WSS profiles for all subjects and the calculated mean blood flow, WSS and peak center velocity.

Figure 6 shows the lumen vessel wall position from one subject during peak systole and end-diastole. The pattern of lateral pulsation in peak systole, as Figure 6 illustrates, was seen in four of the seven volunteers. The cross-sectional lumen vessel area results from end-diastole and peak systole for all subjects are shown in Figure 7. The mean ± SD vessel area pulsation was 28.0 ± 3.1%.

## Discussion

The present report introduces a new method for the noninvasive determination of blood flow, cross-sectional lumen vessel area and WSS based on circumferential subpixel edge detection. The method applies fundamental hemodynamic theory using a 3DP model for fitting of blood velocity data in the boundary layer. The data were obtained by standard MRI velocity acquisition techniques.

Circumferential edge detection was made continuously throughout the heart cycle by an automated procedure solving the 3DP model for u = 0 in the coordinate system, and WSS was calculated by differentiating for u = 0.

The three assumptions used—1) blood velocity at the vessel wall is zero; 2) parabolic velocity distribution is close to the vessel wall; and 3) arterial lumina are circular in shape—have never been definitely proved. These assumptions are a basic concept of general hemodynamic theory (11,12) and have been used in several studies (6,7,11,17–20). Our data did not contradict these assumptions (Fig. 2a, 3 and 4, Table 3).

The large number of data points (n = 69) and selective inclusion of data points from pixels within the blood stream only (avoiding edge pixels [2,6,7]) gave high precision of the fit with a low subpixel magnitude of standard errors (Fig. 4 and 6, Table 3).

### Method evaluation

The in vitro part of this series with Newtonian fluid in noncompliant conduits at steady, laminar and fully developed flow showed that it was possible to determine the internal area of tubes with extreme factory specifications (8.00 ± 0.01-mm internal diameter), with a mean area error of 0.6% of the theoretical value. For both the low flow (diameter ±95% confidence limits 7.961 ± 0.068) and high flow (8.010 ± 0.048) used, the confidence interval included the factory-specified lumen dimensions. Mean blood flow was determined with a 6.2% error, and the calculated mean values for WSS had an accumulated error of 7.3% (Table 1).

To our knowledge WSS derivation and area construction from primary velocity data using a paraboloid fit have not previously been reported, but the low magnitude of errors and standard errors (Table 1) for the in vitro studies substantiate the feasibility of the 3DP method and suggest a potential for in vivo applications.

Blood flow was blunt in peak systole (Fig. 3) but became paraboloid during diastole, and the profiles were slightly asymmetric. The ideal geometric model assumptions were not fulfilled, but the statistical analysis (with subpixel order of magnitude of standard errors [Table 3]) showed that the 3DP method, under the present conditions, accurately described the velocity profile within the boundary layer.

### Comparison with blood flow and vessel dimension studies

The mean volume blood flow in the CCA was 7.42 ml/s (range 5.58 to 9.58) (Table 2), which corresponds to previous determinations with MRI phase contrast methods of (mean ± SD) 7.6 ± 0.9 ml/s (21) and 6.5 ± 0.9 ml/s (19). The diameter of the CCA as determined by our method was 7.39 mm (range 6.22 to 8.78) during peak systole and 6.47 mm (range 5.40 to 7.67) during end-diastole, giving a diameter pulsation of 14.2%, which is slightly higher than the values found by high resolution ultrasound echo tracking techniques (mean diastolic diameter/pulsation percent): 6.2 mm/9.6% (22); 6.5 mm/9.6% (20); 7.4 mm/9.4% (23). In addition, our technique also allows determination of cross-sectional vessel area, and we found an area change of 25.3% to 32.5%. This introduction of vessel area versus vessel diameter increases the sensitivity and accuracy for dynamic determination of vessel dilation/contraction during the heart cycle because radius enters with the second power in the calculation of area (Table 2).

### Comparison with MRI WSS studies

MRI phase contrast velocity mapping techniques are well suited for in vivo estimation of WSS, and two studies have measured WSS in the abdominal aorta using this technique (6,7). To our knowledge, MRI measurements of WSS in the CCA or arteries of similar size have thus far not been presented. We previously introduced (6) an MRI method based on estimating the blood velocity in the edge pixel (partial volume pixel) and in the adjacent pixel within the blood stream. Oshinski et al. (7) proposed a similar method. These methods are not optimal because of the use of partial volume (edge) pixels (11) and linear curve fitting (11). They are both sensitive to noise because only a few pixels are used. The 3DP method overcomes these limitations.

### Comparison with ultrasound WSS studies

New ultrasound methods have been developed for estimation of WSS in the CCA. By assuming a parabolic velocity profile and applying the Poiseuille concept, Gnasso et al. (13) estimated WSS from one central velocity and vessel diameter measurements. When our data are analyzed in a similar way (Table 2), the peak systolic WSS is underestimated by 50% for peak systole, and the mean WSS is underestimated by 22%. This divergence is due to the assumption of a parabolic profile, which our data show to be not valid. That our absolute 3DP data and the ultrasound data from Gnasso et al. (13) (peak WSS 2.95 N/m^{2}, mean WSS 1.21 N/m^{2}) have approximately the same WSS values is partly due to the unusually high central velocities (97.1 cm/s) found by Gnasso et al. compared with other studies (17–19).

Hoecks et al. (20,24) used a newly developed high resolution ultrasound system and also found a blunt profile in systole. They calculated that the mean WSS was underestimated by 86% if a fully developed (paraboloid) velocity profile was assumed (20). For these reasons we suggest that WSS in the CCA cannot be quantified correctly using the assumption of fully developed parabolic flow.

### Special features of the 3DP method

Apart from the estimation of blood flow, WSS and vessel cross-sectional area, the 3DP data on the dynamic distensibility of the arterial wall can be used for estimating elasticity, compliance (if combined with pressure measurements) and other biomechanical properties of the vessel wall. Also, previously undescribed variables, such as the in vivo boundary layer thickness (we estimated the peak systolic boundary layer to be ∼1.6 mm) and accurate in-plane vessel center motion, can be derived.

### Study limitations

Some limitations of the method should be discussed: 1) The assumptions of symmetry of the flow profile and circularity of the vessel—although acceptable in the present study—are not applicable to most arterial vessel segments because of vessel curvature and entrance effects and subsequently skewed velocity profiles with varying boundary layer thicknesses. This problem could be overcome in the future by dividing the whole circumference of a vessel into small sectors evenly spaced around the circumference. With this approach it would no longer be necessary to assume a uniform velocity distribution and a uniform circular shape around the entire circumference of the vessel. As an example, the circumference can be divided into eight sectors each containing 45° of the area. The edge for each sector could then be found as a circular fit in each of these eight sectors. We previously presented (25) preliminary results using this “sectored 3DP method,” and the results are encouraging. 2) The boundary layer thickness is known to vary throughout the heart cycle (11), and a very thin boundary layer should be expected especially in early systole. A quantitative approach for solving this problem would be to estimate the boundary layer thickness for each heart phase and adjust the fit layer accordingly. For end-diastole this could actually mean an increase in the thickness of the fit layer (Fig. 4). 3) The method is somewhat limited by the in-plane resolution and quality of the MRI velocity data. For the present study, the use of standard MRI hardware and software was sufficient, but for smaller vessels and thinner boundary layers and for using the proposed “sectored 3DP method,” higher resolution would improve the accuracy of the method by increasing the number of data points for the 3DP fit. This could be obtained through the use of high resolution techniques, such as Fid Acquired Echo (FAcE) (3), and by using scanner hardware with improved gradient performance, as well as dedicated surface coils (26).

When these current limitations are overcome by the suggested improvements, and a fully automatic (27) and user-independent implementation of the technique is made, the 3DP method can be applied to high resolution MRI velocity data from arteries throughout the vascular system independent of assumptions of velocity profile symmetry and circularity of vessels. Combined with high resolution MRI imaging of plaque morphology (28), this would make a very powerful noninvasive tool with sufficient accuracy for studying intraindividual and temporal development of vascular disease and vascular responses to interventions.

### Conclusions

We described a new noninvasive method for accurate estimation of blood flow, cross-sectional lumen vessel area and WSS. To our knowledge, this is the first method to obtain robust subpixel resolution for cross-sectional lumen vessel area and WSS estimations. In vitro results and statistical analysis demonstrated the feasibility of the method, and the first in vivo implementation using standard MRI hardware and software showed results comparable to published data as well as a low subpixel magnitude of standard errors.

## Acknowledgements

The help of Steen Demming for implementation of software is gratefully acknowledged. Erling Falk is thanked for review of the manuscript.

## Footnotes

☆ The study was supported by the Kirsten Anthonius Foundation, the Danish Heart Foundation, the Karen Elise Jensen Foundation, the Danish Research Academy and the Desirée and Niels Yde Foundation, Copenhagen, Denmark.

- Abbreviations
- MRI
- magnetic resonance imaging
- 3DP
- three-dimensional paraboloid
- WSS
- wall shear stress
- Re
- Reynold’s number
- CCA
- common carotid artery

- Received November 11, 1997.
- Revision received March 26, 1998.
- Accepted April 9, 1998.

- American College of Cardiology

## References

- 1.↵
- 2.↵
- 3.↵
- 4.↵
- Walker P.G.,
- Oyre S.,
- Pedersen E.M.,
- Houlind K.,
- Guenet F.S.,
- Yoganathan A.P.

- 5.
- Kim W.Y.,
- Walker P.G.,
- Pedersen E.M.,
- et al.

- 6.↵
- 7.↵
- 8.↵
- Mohiaddin R.H.,
- Underwood S.R.,
- Bogren H.G.,
- et al.

- 9.
- 10.↵
- 11.↵
- 12.
- Nichols W.W.,
- O’Rourke M.F.

- 13.↵
- Gnasso A.,
- Carallo C.,
- Irace C.,
- et al.

- 14.
- Gnasso A.,
- Irace C.,
- Carallo C.,
- et al.

- 15.↵
- Press W.H.,
- Flannery B.P.,
- Teukolsky S.A.,
- Vetterling W.T.

- 16.↵
- 17.↵
- Enzmann D.R.,
- Ross M.R.,
- Marks M.P.,
- Pelc N.J.

- 18.
- 19.↵
- 20.↵Hoeks AP, Samijo SK, Brands PJ, Reneman RS. Assessment of wall shear rate in humans: an ultrasound study. J Vasc Invest 1995;1:3:108–17.
- 21.↵
- 22.↵
- Reneman R.S.,
- van Merode T.,
- Hick P.,
- Hoeks A.P.

- 23.↵
- 24.
- 25.↵Oyre S, Ringgaard S, Kozerke S, et al. Quantitation of circumferential subpixel vessel wall position and wall shear stress by multiple sectored three-dimensional paraboloid modeling of velocity encoded cine MR [abstract]. Proc Annu Meet Int Soc Magn Reson Med 1998;1:44.
- 26.↵
- 27.↵Kozerke S, Pedersen EM, Scheidegger MB, Boesiger P. Fast flow quantification using a combination of single breath-hold measurements and automated segmentation [abstract]. Proc Annu Meet Int Soc Magn Reson Med 1997;3:1875.
- 28.↵
- Toussaint J.F.,
- LaMuraglia G.M.,
- Southern J.F.,
- Fuster V.,
- Kantor H.L.