Mathematical Models for Blood Flow Quantification in Dialysis Access Using Angiography: A Comparative Study

Abstract

Blood flow rate in dialysis (vascular) access is the key parameter to examine patency and to evaluate the outcomes of various endovascular interve7ntions. While angiography is extensively used for dialysis access–salvage procedures, to date, there is no image-based blood flow measurement application commercially available in the angiography suite. We aim to calculate the blood flow rate in the dialysis access based on cine-angiographic and fluoroscopic image sequences. In this study, we discuss image-based methods to quantify access blood flow in a flow phantom model. Digital subtraction angiography (DSA) and fluoroscopy were used to acquire images at various sampling rates (DSA—3 and 6 frames/s, fluoroscopy—4 and 10 pulses/s). Flow rates were computed based on two bolus tracking algorithms, peak-to-peak and cross-correlation, and modeled with three curve-fitting functions, gamma variate, lagged normal, and polynomial, to correct errors with transit time measurement. Dye propagation distance and the cross-sectional area were calculated by analyzing the contrast enhancement in the vessel. The calculated flow rates were correlated versus an in-line flow sensor measurement. The cross-correlation algorithm with gamma-variate curve fitting had the best accuracy and least variability in both imaging modes. The absolute percent error (mean ± SEM) of flow quantification in the DSA mode at 6 frames/s was 21.4 ± 1.9%, and in the fluoroscopic mode at 10 pulses/s was 37.4 ± 3.6%. The radiation dose varied linearly with the sampling rate in both imaging modes and was substantially low to invoke any tissue reactions or stochastic effects. The cross-correlation algorithm and gamma-variate curve fitting for DSA acquisition at 6 frames/s had the best correlation with the flow sensor measurements. These findings will be helpful to develop a software-based vascular access flow measurement tool for the angiography suite and to optimize the imaging protocol amenable for computational flow applications.

1. Introduction

The ability to measure blood flow using digital images is of significant clinical interest that permits management, diagnosis, and treatment planning of vascular disease [1,2,3]. Autogenous fistula and prosthetic grafts are the most common vascular access types in end-stage renal disease (ESRD) patients that provide “access” to systemic circulation for high-efficiency dialysis [4,5,6]. Because of artificial shunting (fistula/graft), the fistula vein registers an abnormal surge in blood flow and invokes a complex pathophysiological process to accommodate the neo-fluidic stress [7,8]. This results in vascular remodeling (intimal hyperplasia) that, over time, leads to narrowing of the lumen and predisposes to stenosis and thrombosis [9,10,11,12]. Corrective angioplasty (percutaneous transluminal angioplasty) restores the hemodynamic function [13]. As per the National Kidney Foundation Kidney Disease Outcomes Quality Initiative (NKF/KDOQI), the minimum flow threshold to ensure patent access is 400–500 mL/min for arteriovenous fistula and 600 mL/min for arteriovenous grafts [14]. Based on these guidelines, the measurement of post-procedural blood flow can provide a reliable assessment of the radio-surgical endpoint.

The estimation of dialysis access blood flow is currently possible in all imaging modalities, including magnetic resonance (MR) [15,16], computed tomography (CT) [17], or ultrasound imaging [18,19]. The high costs of MR and CT scanning and the related logistic factors prevent its routine use for clinical examination and insurance reimbursement purposes. While the ultrasound technology is relatively inexpensive and does not involve ionizing radiation, multiple factors contribute to the poor estimation accuracy with ultrasound imaging. The error stems from the inaccuracy of the Doppler angle measurement of beam alignment with the vessel’s interface [20,21]. A deviation of > 20° in Doppler angle measurement leads to over 50% error in flow estimation [18]. The poor signal-to-noise ratio limits the ability to accurately measure the cross-section—a factor critical for flow estimation [22]. The attenuation difference between the blood and the tissue is also a concern that causes non-uniform insonation and scattering of the beam frequency [23,24]. Further, the accuracy and interpretation of the ultrasound data depend on the probe type, operator’s skill set, and experience, making the examination subjective [25,26]. Conversely, C-arm fluoroscopy is extensively used for intra-procedural guidance and imaging during dialysis access intervention and has the potential to quantify flow. However, to date, commercial application to measure access flow using fistulagrams is not yet available. The convenience of having an integrated flow measurement application in the angiography suite enables a cost-effective solution and provides an objective assessment of various endovascular procedures.

There have been various efforts to model the circulatory system to study the vascular mechanics and pathophysiological conditions in cardiovascular disease [27,28,29,30]. Polanczyk and colleagues introduced a novel artificial circulatory model (ACM) to simulate hemodynamics under various physiological flow conditions and measured the mechanical response of the vascular tissues [29]. The ACM- an ex vivo computer-controlled bioengineered reactor, was designed to mimic the human circulatory system to study vessel’s wall compliance in real-time. To effect, an “artificial heart” was presented to independently vary the ejection pressure, ejection volume, and frequency of pulsation, and a novel vision acquisition system was conceived to record the 3D-geometrical changes of the vessel in response to different flow stimuli. Vascular biomechanics were simulated under various physiological flow conditions and the corresponding wall deformation responses were validated with the medical data reconstructed with 2D-speckle tracking technique that showed good agreement [29]. In another study from the same investigator, a novel vision-based system was proposed to spatially assess abdominal aortic aneurysm and its wall deformation in a water environment [30]. The non-invasive vision-based system (NIVBS) consisted of a set of nine cameras oriented in circular geometry and custom (self-made) computer vision algorithms were applied to interrogate aneurysm sac wall deformation. Aortic biomechanics were reconstructed in a water environment to simulate physiological pressures from surrounding organs, and the aneurysm diameter and wall deformation were calculated. These experimental data were correlated with AngioCT and 2D speckle–tracking echocardiography that showed good correlation [30]. Xaplanteris et al. have proposed a continuous thermodilution-based approach to calculate absolute coronary blood flow (Q) and microvascular resistance (R) in human subjects [31]. The authors used a monorail thermodilution catheter equipped with temperature and pressure sensors to measure the flow and resistance by injecting a saline bolus at room temperature. The quantification of microcirculation and resistance allows a simple and operator-independent assessment and treatment of coronary microcirculatory dysfunction in patients with angina. The absolute coronary blood flow and resistance were calculated based on the following equation [31]:

$$Q=c_p\left[\frac{\left(T_b-T_i\right)}{\left(T_b-T\right)}-1\right]Q_i+Q_i$$

(1)

where, Q_i is the saline infusion rate (mL/min), T_b is the temperature of blood before infusion (≈37 °C), T_i is the temperature after infusion at catheter exit, T is the homogeneous temperature of the blood-saline mixture, and c_p is a constant related to the thermal properties of blood and saline (=1.08). Upon simplification and eliminating negligible quantities, Q can be represented as [31]:

$$Q=1.08\frac{T_i}{T}{Q}_i$$

(2)

R was calculated following Ohm’s law analogy:

$$R=\frac{P_d}{Q}$$

(3)

where, P_d is the pressure at the distal coronary.

Medical image analysis allows extraction of the data that are normally buried in the images and application of relevant mathematical models to extract the physiological parameters of interest. These hemodynamic metrics, such as velocity, pressure/stress, wall diameters, etc., can be correlated to the disease pathology and can foster our understanding of disease initiation/progression and seek rational solutions. Khanmohammadi and colleagues have proposed a multi-step automatic method to estimate blood flow velocity in coronary arteries based on cine X-ray angiographic sequences [32]. X-ray coronary angiography is an invasive routine procedure performed to assess coronary artery disease in patients with angina pectoris. The study recruited 21 patients who did not have obstructive coronary artery disease and coronary flow velocity was measured. The process of velocity calculation involved segmentation of coronary arteries and area calculation, removal of heart motion, polynomial fitting (7th order) on area sequences for start/stop time determination of indicator passage, and velocity calculation based on the principles of indicator propagation. The calculated flow velocity significantly correlated with the transthoracic Doppler recordings [32]. Similarly, intracoronary frequency-domain optical coherence tomography (FD-OCT) was used to examine blood flow rate and velocity in coronary artery stenosis [33]. The measurements were correlated with fractional flow reserve (FFR), which is considered a gold standard for the evaluation of coronary stenosis. OCT uses near-infrared light as opposed to sound waves by intravascular ultrasound (IVUS) and provides high-resolution imaging for lumen measurement. The study evaluated 20 coronary stenoses in 15 patients successively with quantitative coronary angiography (QCA), FFR, and FD-OCT. The authors reported a statistically significant correlation between FD-OCT and FFR calculated flow velocity that suggested the potential value of FD-OCT for the assessment of coronary artery disease against QCA and IVUS imaging modalities [33]. In our study, we are interested to quantify blood flow in dialysis access based on the tools available in the radiological suite. Without additional logistic requirements, the software-based approach, when fully validated, is envisioned to seamlessly integrate into the existing image analysis suite as an installable add-on package.

The purpose of this study is to evaluate image-based flow quantification methods for dialysis access using bolus tracking algorithms and curve-fitting models and to inform the most pertinent method applicable for flow quantification. In our previous study, we quantified blood flow based on the images of digital subtraction angiography (DSA) acquisition and did not employ image pre-processing or curve-fitting models that resulted in poor reproducibility and inaccuracy of measurement [34]. Here, it is hypothesized that the use of curve fitting reduces variabilities in time measurement and results in an enhancement of accuracy and consistency of measurement.

2. Materials and Methods

A radiographic technique to measure blood flow in artificial dialysis access is presented here. Blood flow measurement using angiography is based on tracking the passage of the radiopaque (iodinated) contrast material (CM) at two regions of interest (ROIs) and measuring the corresponding transit time, spatial displacement, and vessel dimension [35]. Using a flow phantom model and a C-arm fluoroscope, hemodynamics from the injection of CM was studied.

2.1. Transit Time Algorithm

The transit time algorithm determines the time taken by the contrast bolus to travel between two user-specified ROIs. A time-dependent intensity map of the transiting bolus is generated at each ROI, also known as time-density curves (TDC), provides the basis for time calculation. Based on various calculation methods (algorithms), the characteristic feature(s) of the first curve is tracked and compared with the second curve [36,37]. Herein, we discuss two major bolus tracking algorithms for transit time estimation. A comprehensive review of other applicable transit time algorithms has been discussed elsewhere [38,39,40].

2.1.1. Peak-to-Peak (PP) Algorithm

The PP algorithm calculates the bolus arrival time based on the time corresponding to the appearance of peak (maximum) intensity at each ROI. The transit time (τ_PP) is then calculated as the difference between the bolus arrival times at each ROI [35].

[t_ROIx (x = 1, 2)]_PP = t(P_ROIx)_max

(4)

$${\tau }_{PP}=\Delta t={\left[t_{ROI2}\right]}_{PP}-{\left[t_{ROI1}\right]}_{PP}$$

(5)

where, t_ROIx is the time corresponding to the appearance of maximum peak intensity (P)_max at ROI_x, and Δt (=τ_PP) is the required bolus transit time.

2.1.2. Cross-Correlation (CC) Algorithm

The CC algorithm calculates the transit time based on a process known as template matching, i.e., the TDC obtained at one ROI is compared with the TDC at another ROI [34,41]. The process of CC requires time-shifting of one curve against the other until a maximum overlap or an optimal match is obtained. The time for which the CC function (φ) is maximized is the required transit time or the initial delay between the two curves.

$${\phi }_{P1P2}\left(\tau \right)=\frac{1}{N}{\displaystyle \sum }_{t=1}^N{P}_1\left(t\right)P_2\left(t-\tau \right)$$

(6)

where P₁ and P₂ are the time-density functions corresponding to ROI1 and ROI2, respectively, N is the number of temporal data samples, and τ = N × Δt; Δt is the unit time interval (time resolution) of the imaging modality.

2.2. Curve-Fitting Models

The calculation of the transit time using the raw signal data is prone to errors from random fluctuations of intensity, quantum noise from image acquisition (image mottle) [42], and recirculation [43]. The curve-fitting function provides a mathematical description of the TDCs based on the statistical properties and assumptions of the underlying distribution. In theory, the curve-fitting function reconstructs a noise-free reference curve amenable for computational purposes by isolating the bolus first pass, which can be used for the estimation of transit time. In this study, we evaluated three curve-fitting models.

2.2.1. Gamma-Variate (GV)

The GV function is one of the most commonly used curve-fitting methods to model the indicator-dilution phenomena in hemodynamic studies [15,19,44], first described by Thomson et al. [45]. It replaces the original intensity curve with a characteristic curve that has a sharp ascending slope, a slowly decaying tail, and a compact spread (narrow curve width) that closely resembles a right-skewed Gaussian curve. The fitting parameters can be optimized to ensure the best fit to represent a particular phenomenon.

Mathematically, the GV function can be written as [45]:

$$C\left(t\right)=K_a{\left(t-AT\right)}^{\alpha }{e}^{-\left(t-AT\right)/\beta }$$

(7)

where, C(t) is the indicator concentration at time, t, K_a is the constant scale factor that determines the level (maximum magnitude) of C(t), AT is the contrast appearance time, and α and β are arbitrary parameters that determine shape (rate of decay) of the curve.

2.2.2. Lagged Normal

The lagged normal distribution was first described by Bassingthwaighte et al. to model the dispersion of an indicator in the flowing fluid [46]. It was derived to overcome Poiseuille’s flow model (parabolic flow) for estimating blood flow in the vasculature, which is often skewed and turbulent. It was shown that the model can fit a wide variety of dilution curves. Lagged normal function describes the indicator-dilution phenomenon based on the compartmental model of tissue perfusion that combines the normal distribution (Gaussian function) with an exponential function to represent the trailing edge [47]. In particular, the Gaussian function provides the basis for random dispersion with turbulence, and the exponential function models the mixing that occurs in the chambers [48]. Lagged normal distribution was used to model the blood flow from a large vessel into a microvascular bed [19].

In the simplest form, we can represent lagged normal distribution by the convolution of two functions: a Gaussian function and an exponential function [19,46,48].

$$C\left(t\right)=\left[f\ast g\right]\left(t\right)={{\displaystyle \int }}_{-\infty }^{\infty }f\left(\tau \right)g\left(t-\tau \right)d\tau ={{\displaystyle \int }}_{-\infty }^{t}f\left(\tau \right)g\left(t-\tau \right)d\tau$$

(8)

(Since, g (t − τ) = 0 for t < τ, the upper limit of the integral is t)

Where

$$f\left(t\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-{\left(t-\mu \right)}^2/\left(2{\sigma }^2\right)},-\infty <t<\infty$$

(9)

$$g\left(t\right)=\lambda e^{-\lambda t}, t\ge 0$$

(10)

where µ and σ² represent the mean transit time and variance of the Gaussian function. The parameter λ is the rate constant of the exponential function.

Solving for C(t) resolves into:

$$C\left(t\right)=\frac{A}{2}K\left[1+\mathrm{erf}\left(L\right)\right]$$

(11)

where A is the area under the curve, and erf(.) indicates the error function, defined as:

$$erf\left(x\right)=\frac{2}{\sqrt{\pi }}{{\displaystyle \int }}_0^x{e}^{-t^2}dt$$

(12)

The parameters K and L are given as:

$$K=\lambda e^{\left(-\lambda t+\lambda \mu +\frac{1}{2}{\lambda }^2{\sigma }^2\right)}$$

(13)

$$L=\frac{t-\mu -\lambda {\sigma }^2}{\sqrt{\left(2{\sigma }^2\right)}}$$

(14)

2.2.3. Polynomial (6th Degree)

The lower degree polynomial fails to fit the raw data adequately to resemble a typical TDC. Conversely, higher-order polynomials fit the data quite well along with noise, which is undesirable. Therefore, the process of fitting and selection (optimization) of polynomial order requires a cautious approach to balance the accuracy of fit against noise inclusion.

The sixth-order polynomial, used in this study, can be written as:

$$C\left(t\right)=a_1+a_{2}t+a_3{t}^2+a_4{t}^3+a_5{t}^4+a_6{t}^5+a_7{t}^6$$

(15)

$$C\left(t\right)={\displaystyle \sum }_{i=1}^7{a}_i{t}^{i-1}$$

(16)

where C(t) is the contrast concentration at any time, t, and a_i is the fitting constant to be determined.

2.3. Optimization of Curve-Fitting Parameters

Initially, an approximated (coarse) value of the unknown parameters was obtained by plugging random numbers. Based on the direction of the reduction (minimization) of the cost function (root mean square error, RMSE), the parameters were further optimized in fine increments. As an example, we have provided the process of optimization for the GV curve-fit function (Figure 1). The results of other curve-fitting models are depicted. (Note: the polynomial fitting was executed using MATLAB’s inbuilt “polyfit” function).

2.4. Distance and Cross-Sectional Area Calculation

The edge detection method was specifically developed to calculate the spatial distance traveled by the contrast bolus and the vessel’s cross-sectional area [34,35]. Because the injection of the CM enhances the vessel (contrast opacification), edges can be determined by discriminating the intensity interior of the vessel with the image background. The algorithm for edge detection is illustrated (Figure 2). For every row, assume each box in the image to be of unit pixel length of given intensity (white pixel- low intensity, black pixel- high intensity, gray pixel- in between intensity (sense of pixel intensity was inverted to guide the eye)). The assignment of a row begins with the coordinates of ROI1 and ends with the coordinates of ROI2. Each row within the ROI produces exactly one distance-intensity plot. In this manner, multiple distance-intensity curves were plotted along the vessel’s length (bolus travel path) at discrete intervals (25-pixel separation). At the vessel interface (either edge), there was a sharp change in the contrast bolus intensity. The initial cut-off intensity to declare edge transition was set at 5% of the maximum opacification and fine-tuned until an accurate edge identification could be reliably made (settled at ~12%). After the edges were detected, the centerline was calculated as the mean distance between the two extracted edges. The contrast traversal distance was calculated based on the pixel separation between the two ROIs along the extracted centerline. Because the image mottle was more prevalent on fluoroscopic images, the extracted edges and centerline were relatively noisier over DSA acquired images.

The cross-sectional area was calculated assuming a circular geometry. The access diameter was determined from the curve width of the distance-intensity plot, generated previously (edge detection). To accommodate for changing vessel’s diameter often encountered in clinics, the mean cross-sectional area approach was adopted for flow estimation, i.e., the average of multiple cross-sectional areas sampled at discrete intervals (25-pixel separation) between the ROIs was used.

2.5. Flow Phantom Set Up

A water-based benchtop model for simulating dialysis access was constructed (Figure 3). It comprised a water tank, a peristaltic pump (EW-07528-10 & EW-77200-50, Cole-Parmer Instrument Company, Vernon Hills, IL, USA), power injector (Mark V Provis, Medrad Inc., Warrendale, PA, USA), in-line flow sensor (ME6PXN, Transonic Systems Inc., Ithaca, NY, USA), and a discard serially connected in an open circuit configuration via a silicone tubing [34,49]. Additionally, three different tubing configurations (straight, angular, and loop) were tested. The entire setup was placed over an angiographic table and imaged in DSA and fluoroscopy filming modes using a C-arm fluoroscope.

The DICOM images obtained from the above imaging experiments were fed offline into an ordinary personal computer equipped with a GUI-based access flow computational model constructed entirely in MATLAB. It required the users to identify two ROIs in the bolus flight path and had the option to use various curve-fitting models and interpolation functions for flow calculation.

2.6. Image Acquisition Protocol

The imaging protocol was developed to resemble a realistic clinical imaging environment [34,49]. DSA images were acquired at 3, and 6 F/s; fluoroscopic images were acquired at 4, and 10 pulses/s (P/s). The total duration of imaging was set at 8 s, and CM was power injected via a sidearm 6 French introducer sheath. The contrast injection was synchronized with the imaging function in the C-arm fluoroscope and injected at 3.33 cc/s for 3 s (10 mL over 3 s) with an initial delay of 1 s. The imaging duration and initial injection delay were optimized based on computational experience that required a pre-contrast phase along with bolus wash-in and wash-out segments for enhanced accuracy.

3. Results

The overall results obtained with each algorithm with and without different curve-fitting models were plotted to examine a general trend of accuracy and reproducibility for each imaging mode (Figure 4). We quantified the accuracy in terms of absolute percent quantification error, i.e., absolute percent deviation from in-line flow sensor measurement (ground truth) for each observation. For both the imaging modes, the CC algorithm had the highest accuracy (least quantification error) and least deviation (high precision). Among the different curve-fitting models, the GV curve fit had the best accuracy and precision of measurement. For DSA acquisition, the lowest quantification error (mean ± SEM) was obtained with CC + GV curve fit (6 F/s), and the highest quantification error of 52.9 ± 4.3% was obtained with the PP algorithm without curve fit (4 F/s). Since the CC algorithm had the best overall accuracy, the effect with/without curve fit on accuracy at different sampling rates was evaluated (Figure 4c). The quantification error with GV curve fit at 4 F/s was 27.6 ± 2.7% and at 6 F/s was 21.4 ± 1.9% (highest accuracy). Similarly, for fluoroscopic imaging, the least quantification error (mean ± SEM) was obtained with the CC algorithm + GV curve fit (10 P/s), and the highest quantification error of 64.1 ± 6.5% was obtained with PP algorithm without curve fit (4 P/s). The variation in quantification accuracy with the CC algorithm at different pulse rates as a function of curve fit or none (no curve fit) was plotted (Figure 4f). The quantification error with the GV curve fit at 4 P/s was 44.4 ± 4.3% and at 10 P/s was 37.4 ± 3.6% (highest accuracy). The quantification error of the PP algorithm with different curve-fitting models for both imaging modalities (DSA acquisition and fluoroscopy) is illustrated (Figure S1). The quantification error with the GV curve fit at 6 F/s was 26.3 ± 2.9% and at 10 P/s was 42.2 ± 4.3%. The variation in quantification accuracy with the CC algorithm as a function of flow rate for DSA acquisition is depicted (Figure 5). The results are shown for the CC algorithm with or without GV curve fit because of the consistency of measurement, as stated previously (for PP algorithm, refer to Figure S2). Data show DSA acquisition at 6 F/s with CC algorithm + GV curve fit had the least quantification error when the flow rates were varied. These results affirmed that the CC algorithm with GV curve fit had the highest accuracy and precision of measurement and is best suited for flow computational studies.

The bias, limits of agreement (LOA), and 95% confidence interval (CI) of the mean quantification error for the DSA acquisition and fluoroscopic imaging were assessed (

Table 1. Summarized key statistical parameters for PP and CC algorithms with/without curve fitting.

Imaging Mode, Sampling Rate	Algorithm	Without Curve Fit			With GV Curve Fit			p Value †
		Bias [Lower Limit, Upper Limit] (mL/min)	95% CI |MQE| (%)	Corr.* (r)	Bias [Lower Limit, Upper Limit] (mL/min)	95% CI |MQE| (%)	Corr.* (r)
DSA acquisition, 6 F/s	PP	−176 [−842, 490]	[29, 43]	0.574	−246 [−676, 184]	[21, 32]	0.773	0.175 ‡
	CC	−234 [−691, 223]	[25, 35]	0.731	−145 [−484, 195]	[18, 25]	0.865	0.017 §
Fluoroscopic, 10 P/s	PP	−374 [−1098, 351]	[42, 67]	0.545	−339 [−891, 213]	[34, 51]	0.515	0.027 ‡
	CC	−205 [−930, 519]	[43, 59]	0.534	−198 [−722, 327]	[30, 44]	0.576	0.031 §

* Pearson’s correlation coefficient (r) of calculated flow rate vs. in-line flow sensor measurement. ^† Statistical test of significance (Student’s t test) between algorithms (peak-to-peak vs. cross-correlation; ^‡ without curve fitting, ^§ with gamma-variate curve fitting, respectively).

). Bias and LOA were calculated based on Bland-Altman’s theory [50]. Bias indicates the systemic tendency of the method to under- or overestimate the actual measurement and calculated as the mean of the measured differences with in-line flow sensor measurement. LOA indicates agreement of the calculated value with the measured (actual) value and depicted by a bandwidth (bias ± 1.96 × SD of differences), where 95% of the differences are expected to lie (the narrower the bandwidth, the more precise is the measurement). Since the CC algorithm had the least quantification error, we have shown BA analyses for this method (DSA acquisition: Figure 6, Fluoroscopy: Figure S3). The CC algorithm with GV curve-fitting had the best agreement with the measured values. The 95% CI for the bias and LOA bandwidth was also the least indicating reproducibility (precision/ repeatability) of the measurements. The CC algorithm with GV fit at 6 F/s exhibited a good correlation with the actual flow (r = 0.898). We made similar observations for fluoroscopic imaging, except that the range was wider, which indicated large variabilities of measurements over DSA acquisition. The correlation (r) between measured flow rates obtained with CC algorithm + GV fit at 10 P/s and in-line flow sensor was 0.334.

Since the use of C-arm angiography involves ionizing radiation, the radiation burden at different sampling rates was evaluated for DSA and fluoroscopic examinations (Figure S4). In clinics, radiation doses are usually stated in terms of air kerma (AK) and air kerma area product (AKAP, or dose area product (DAP)) since these quantities are readily available in the modern fluoroscopic machines (angiographic workstation) and collected at the end of the procedure for record-keeping purpose. These parameters can be used to derive the peak skin dose (PSD) and effective dose (ED), which are more sensitive and specific quantities for measuring the effects of radiation, i.e., tissue reactions (deterministic effect (PSD)) and stochastic effects (radiation-induced cancer (ED)). It was found that the radiation dose from these procedures was substantially low to invoke any significant detrimental effects. Radiation dose varied linearly and increased proportionately with the sampling rate.

4. Discussion

This study aimed to identify key algorithm(s) and curve-fitting model(s) that can provide an accurate and consistent estimation of the access flow should a software application be developed. We tested two transit time algorithms and three curve-fitting models to find an optimum method for quantifying dialysis access flow. It was found that the CC algorithm with GV curve fit, compared to other methods, had the most accurate and consistent description of the access flow (DSA—6 F/s, Fluoroscopic—10 P/s). The CC method uses the area under curve approach to calculate the optimum time delay that makes it comparatively immune to random fluctuations of the intensities [34,41,49].

During the radiological intervention of the dialysis access, a concentrated bolus of CM is injected into the bloodstream that enhances the path of the bolus flow [51,52,53]. While the visibility of the vascular network is fundamental for clinical diagnosis, it also contains intrinsic hemodynamic information that can be exploited for measuring physiological flow [34,35]. The estimation of bolus arrival time using raw data is challenging, as it contains noisy signals and is devoid of any characteristic shape [49]. Noisy components appear on the dilution curves because of the non-homogenous composition of the blood-contrast mixture and technological factors associated with imaging [54]. Since bolus arrival time depends on the shape of the enhancement curve, fluctuations of intensities change the curve contour and exert a negative impact on computational accuracy. Therefore, a uniformly enhanced curve profile with minimum intensity fluctuations is expected for quantification, which can be accomplished through the use of curve fitting. The GV curve-fit function had the best results among all other curve-fitting methods. It provided a close similarity to the actual physiological flow with a relatively simpler optimization process. The lagged normal distribution could also model the dilution curves and get rid of unwanted random intensity outliers; however, the accuracy and repeatability were lower than that of the GV curve-fitting approach. Conversely, polynomial fitting did not accurately model the TDC and more likely followed the noise inherent in the dilution curves that resulted in poor computational accuracy.

During clinical imaging, the sampling rate and imaging duration are usually kept low to avoid a substantial amount of radiation dose to the patient and the clinical staffs in the operating room [55,56,57]. The increase in sampling rate affects the computational accuracy and the radiation dose as it relates to the number of samples available for computational purposes and the X-ray beam on time, respectively. At a low sampling rate, observations were obtained at infrequent intervals that amounted to a low number of data samples and poor depiction of the bolus passage. At a high sampling rate, large amounts of data (better image quality) were rapidly generated, albeit at the cost of increased radiation burden. Thus, a balance between accuracy and acceptable radiation dose needs to be ascertained for a clinically realistic environment. The radiation burden increased linearly with the increase of frame rate or pulse rate but was substantially low to invoke any radiation-induced injury. The sampling rate of 3 F/s (DSA) or 4 P/s (fluoroscopy) had the least accuracy. With the increase in sampling rate to 6 F/s or 10 P/s, computational accuracy was improved. The quantification accuracy and precision with the DSA acquisition were superior to that of the fluoroscopic imaging. The evaluation of fluoroscopic imaging was motivated by its extremely low radiation footprint. However, because of the use of a low radiation X-ray beam, the images were noisier that affected the extraction of the cross-section area and transit time. The accurate estimation of the cross-sectional area is of paramount significance as it impacts the flow calculation to the second power of the radius (∝ R²). Likewise, the conversion of flow rate units (mL/s to mL/min) involves a factor of 60 that rapidly overestimates the calculated flow rate, including with a small degree of transit time calculation inaccuracy. Therefore, for clinical purposes, DSA acquisition at 6 F/s may be preferred. Waechter et al. used 3D rotational angiography to model the flow and to determine hemodynamic parameters based on the spatial and temporal progression of the contrast concentration in a phantom model [3]. The error rate was 5–10% of actual flow, and the images were obtained at 30 F/s. Similarly, Pereira and colleagues used 3D rotational angiography to study intra-aneurysm hemodynamics induced by flow diverter stents that led to progressive thrombosis and vascular remodeling [58]. The study was conducted on 24 patients with variable flow rates. Images were acquired at 60 F/s, and the correlation of flow amplitude with the prediction of aneurysm thrombosis was 88% (sensitivity) and 73% (specificity). A cross-correlation algorithm was used for flow quantification by Seifalian et al. [59]. Images were acquired at 25 F/s; the error of quantification was 20–50% for flow velocity ranging from 147 to 733 mm/s. Against this background, this study obtained an accuracy of 21.4 ± 1.9% with DSA acquisition at a low sampling rate of 6 F/s, which is a realistic frame rate for clinical applications.

The measurement of access flow involves a complex interplay of multiple factors that includes optimizing the parameters of injection and imaging duration and managing the radiation dose in a manner that maximizes the diagnostic value with minimum risk for radiation-induced injury [43,54,60,61]. Since the determination of blood flow is entirely dependent on the analysis of radiological exams, it is necessary to have good quality images with adequate temporal information to enable accurate quantification [62,63]. With these facts in mind, we developed a custom injection and imaging protocol for use with this study. The objective of the customized protocol was to record the pattern of physiological flow during image acquisition, possibly with the lowest radiation burden. We emphasize the use of a power injector over the handheld injection method for administering CM as human factors such as the rate and volume of injection can affect flow dynamics. The use of a power injector produces highly uniform enhancement curves without operator dependency and delivers precision in time calculation. In addition, the concentration of the CM directly affects the quality of enhancement. In this study, we used CM at full strength (100%) and did not experiment with various dilutions since the volume of administration was considerably low to invoke any adverse reactions to CM non-allergic patients [64]. However, while the use of diluted or concentrated CM media may have no substantial impact on the analysis of DSA images, it can affect the analysis of fluoroscopic images, as the image quality is extensively degraded (poor signal-to-noise ratio).

We evaluated contrast enhancement curves to calculate the centerline and distance traveled between two regions of interest. We first applied this concept to a simple tube configuration followed by angled and loop configurations and could extract the centerline and distance traveled by the contrast bolus with reasonable accuracy.

This study has few limitations. First, we used water in lieu of blood and it may not necessarily correlate with the in vivo fluid dynamics. For instance, the hematocrit influence on blood flow was not studied. However, while the hematocrit does affect viscosity, the system was conceived as an internally controlled model with a rigorously validated gold standard (ultrasonic flow sensor) for the measurement of flow. Because we measure the same material with angiographic imaging and ultrasonic flow sensor, the composition of the flowing fluid may not exert a significant impact. Further hematocrit effects on flow mostly relate to capillary flow and the physiological effects of altering the cellular composition of the blood on the apparent viscosity of blood (Fahraeus effect) [65]. The fistula is a clinically constructed conduit in the peripheral circulatory system to facilitate dialysis with vessel’s dimension ranging up to 1–2 cm in diameter and registers a flow rate of 300–1000 mL/min (patent access). Due to the high bulk flow rate in the dialysis fistula, the homogenous blood-contrast mixture is rapidly dispersed in circulation [43]. For more realistic blood flow simulation, future studies need to examine blood analogs such as water-glycerol or xanthan gum-water-glycerol mixtures as potential test fluids [29,66]. Second, we chose a simple geometric construction of the vessel and applied circular cross-sectional model for area calculation. These assumptions were made to simplify the study and to avoid the use of 3D rotational angiography to acquire the third-dimensional parameter necessary for the reconstruction of a 3D model. The assumption of a circular cross-section made here may not necessarily hold for stenotic vessels, which are essentially malformed. Since the change in vessel geometry affects the flow regimen (oscillatory/turbulent flow), future study needs to investigate the use of a noncircular cross-sectional model, its impact on dispersive flow, and the technical aspects associated with image acquisition and analysis.

5. Conclusions

A radiological concept to estimate flow based on the principles of indicator propagation under the influence of bulk flow was presented. The measurement of blood flow provides the functional significance of stenosis and elucidates the severity of the pathology that permits immediate remedies [67,68,69]. The paper introduced and compared algorithms based on their respective accuracy and reproducibility. The estimation of flow rate using the DSA and/or fluoroscopic images is a clinically relevant problem; therefore, the identification of reliable algorithms and methods in this direction is of particular significance. The CC algorithm with GV curve fit had the best overall accuracy for time measurement. For optimal accuracy, DSA acquisition at 6 F/s needs to be used. This study has provided a necessary groundwork to conduct a small-scale clinical study for feasibility analysis based on the approaches and parameters outlined and may lead to the development of a novel imaging product for vascular access application.