Comparison of various zero-dimensional models to analyze the brachial artery blood flow of an arterio-venous fistula during hemodialysis procedure

One of the most recommended vascular accesses for hemodialysis is an autogenous arterio-venous fistula (AVF), which is a surgically created conduit between an artery and a vein. AVFs can be compromised by stenosis-induced arterial or venous thrombosis. An effective and efficient mathematical fluid mechani-cal model can provide insight into the mechanisms for stenosis and thrombosis. In an earlier paper, a laminar entry flow model (Shah) was applied to a fifth order lumped parameter zero dimensional (0D) model to predict the inlet flow rate in the brachial artery of an AVF. ODE45, a built in MATLAB ordinary differential equation (ODE) solver, was used to produce the solution. The Shah model showed closer correlation with a computational fluid dynamics model, performed in COMSOL, than did the Poiseuille flow model. In this study, an outlet flow rate in the brachial artery of an AVF is generated with third, fourth, and fifth order lumped parameter zero dimensional models with laminar and turbulent flow conditions within different arterial segments. Butcher’s sixth order Runge-Kutta (RK) method was used as an ODE solver and the results were compared with the MATLAB ODE solvers (ODE45 and ODE113). The lumped parameter model solved with Butcher’s sixth order RK can help to understand the AVF performance in vitro and guide the application of an appropriate lumped parameter model for blood flow prediction.

can help to understand the AVF performance in vitro and guide the application of an appropriate lumped parameter model for blood flow prediction.

arterio-venous fistula (AVF), blood flow modeling, hemodialysis, lumped parameter INTRODUCTION
Hemodialysis is a life-saving therapy for patients suffering from the end-stage renal disease (ESRD). The total number of ESRD cases in USA in 2016 was over 700,000, of which 87.3% of the incident individuals began therapy with hemodialysis. 1 The incidence of ESRD in Nepal is more than 2900 cases (100 per million population). 2 In hemodialysis, blood is circulated to a dialyzer to remove waste and add nutrients, after which it is returned to the venous side at approximately 300-3000 ml/min. [3][4][5] Permanent vascular access is preferred to carry out the procedure of hemodialysis. Of the many forms of vascular access, a matured autogenous arterio-venous fistula (AVF), which is formed by direct connection of a vein to an adjacent artery, has longest survival rates with fewest complications. 5 About 30%-50% of attempted AVFs fail to mature for hemodialysis. 6 Some studies consider a matured AVF to be one that can deliver a minimum blood flow of 300-400 ml/min in dialysis sessions for prolong time period. 5,7 However, an initially matured AVF is not guaranteed to continue functioning adequately in subsequent dialysis sessions. 5 The exact mechanism of AVF failure remains unresolved, however, some studies claim that local hemodynamic factors such as fluid shear stress and pressure-flow relationships of the blood, could influence AVF maturation by impairing vein dilatation and arterializations and by promoting stenosis and thrombosis. 4,8 Many of the hemodynamic factors of an AVF cannot be practically obtained during a clinical dialysis session. Caroli et al noted that ultrasound signals are often inadequate for blood flow measurement at the AVF site because blood velocity profiles are non-uniform. 7 Therefore, a mathematical model of an AVF may provide insight about AVF maturation to the health care professional, without additional interventions. Other studies have also supported the idea that mathematical model can be useful and may even prolong the life of an AVF and its maturation. [9][10][11][12] A Windkessel (WK) type lumped parameter model can capture various physiological parameters of the brachial artery in an AVF. The model is an electrical circuit-like network where electrical resistance, inductance, and capacitance are analogous to fluid resistance, fluid inertia, and arterial elasticity, respectively. These parameters form a coupled system of ordinary differential equations (ODEs) that mathematically represent blood flow in the brachial artery. The zero dimensional (0D) or lumped parameter model has been used to study the blood flow mechanism for a specific physiological condition. 11 Stergiopulous et al discussed the limitations of second and third order WK models and have introduced inertial term in the fourth order WK model to see the difference between the accuracy of the flow as predicted by different orders of WK models. 13 The shape of the flow waveform (influenced by the pulsatility) is affected by the presence of abnormalities in the arteries. 14,15 A WK model requires the solution of a set of coupled/single ODEs. An analytical solution of these equations has been a challenge, especially when the ODEs are nonlinear. One-step methods, such as the Runge-Kutta (RK) method, are widely used to solve such a set of ODEs numerically. Puelz et al used a second order RK method to solve a third order WK model of1D blood flow in a large vessel network. 16 To understand the hemodynamics of the arterial blood flow in the vascular surgery, Marchandise et al applied the RK method to a patient specific WK model. 17 In our previous study, 4 blood flow in the brachial artery of an AVF was modeled with a fifth order WK model (pulsatile) for different flow conditions (i.e., Poiseuille and Shah). The built-in MATLAB adaptive RK function ODE45 was used. The results were compared with a computational fluid dynamics (CFD) simulation of the brachial artery blood flow performed in the commercially available software COMSOL Multiphysics. Variations in the results from those two studies ranged from 6% to 35%, depending upon the flow conditions. 4 Because arterial stiffness is linked to AVF maturation failure, 4 vessel wall elasticity (pulsatility) is crucial in modeling the brachial artery blood flow.
Dallas et al compared MATLAB ODE solvers to hybrid problems, that is, stiff and non-stiff problems, and found that ODE113 ranked first followed by ODE45. 18 ODE45 is suitable for a linear system of ODEs for pharmacokinetics and pharmacodynamics model. 19 However, Butcher found that higher order RK method were more accurate and time saving when solving the system of ODEs. 20 Blood flow resistance in the brachial artery can be modeled by Poiseuille's law for laminar fully developed flow, and by the Shah model for laminar entry flow. 21 The Shah relationship has closely matched the pressures and flows in brachial artery with an in vitro model of vascular access and with CFD results. 4,12 In the previous study, only laminar flow conditions were considered 4 however, a study by Dixon suggested that the turbulent blood flow might start in the brachial artery, proximal to the anastomosis in an AVF. 5 This claim is further supported by a clinical Doppler ultrasound study, which was conducted by Zamboli et al, where they found blood flow turbulence at the end of the brachial artery in an AVF. 22 In this case irregular flow patterns rather than high Reynolds Number (Re) denoted turbulent flow. Paulson et al stated that there is no strict threshold for turbulence as laminar flow can occur with Reynold Number, Re>2000 and turbulent flow can occur at Re<1000. 23 Therefore, to understand the flow mechanism more properly at the anastomosis of an AVF, Blasius turbulent flow model is used in this study.
For the first part of the analysis, this study uses similar configuration as mentioned in the earlier study. 4 Third and fourth order lumped parameter models remain popular in blood flow modeling. Since, fifth order WK model from the earlier study was validated with the CFD results, in the first part of this study, comparison between the third, fourth, and fifth order lumped parameter model with different ODE solvers is performed to observe how the results differ from each other and to determine whether the ODE solver affects the obtained results. The impact of the ODE solver was not studied in the previous study, therefore, in this study multistep methods such as Adams-Bashfourth-Moulton (ODE113 in MATLAB), ODE45, and a higher order ODE solver, namely Butcher's explicit sixth order RK, were used to solve the ODE, and the results were compared against each other. This study elucidates the influence of the ODE solver on the performance of the ODE (which represents the brachial artery blood flow).In the second part, the results are validated with real clinical data. Blood flow resistant models (i.e., Poiseuille, Shah, and Blasius) are used to model the blood flow at different segments of the brachial artery in an AVF and the result are compared with the clinical data.

METHODS
The geometry of the brachial artery is idealized in Figure 1. The artery length was 20 cm, and the wall thickness was 0.38 mm. For Part-I of the analysis, the inner radius was 2 mm in the unloaded state and ranged from 2.2 to 2.36 mm between diastolic and systolic pressure. The blue region represents the blood flow at the inlet region. The dynamic viscosity and density of the blood were 0.004 kg/(m s) and 1100 kg/m 3 , respectively. 4 For the second part (Part-II), the dynamic viscosity and density of the blood were 0.003 kg/(m s) and 1050 kg/m 3 , respectively. 3,7 These dimensions were extracted from the literature. 7 The fifth order WK model is shown in Figure 2. Figures 3 and 4 represent the fourth and third order WK model, respectively. Kirchhoff's law is applied to all of the models to obtain the required differential equation for the outflow rate. R p is the peripheral resistance and the normal value of R p , before AVF creation (7 × 10 9 Pa s m −3 leads to a flow rate of 117 ml/min, which matches the normal range of brachial artery flow. 4 The lower value of R p , after AVF creation (8.75 × 10 8 Pa s m −3 ) leads to a peak flow rate of approximately 1100 ml/min, representative of an AVF.
The model equations for fifth order WK model ( Figure 2) are: where the currents (i) are analogous to flow rates, and the voltages (v) are analogous to pressures. Here,i 1 is the flow at the inlet of the brachial artery, R 1 is the flow resistance at the inlet of the brachial artery, v in is the pressure at the inlet of the brachial artery, v brach is the pressure at the brachial artery, C brach is the elasticity of the brachial artery, L 1 is the inductance or inertia at the inlet of the brachial artery, i 2 is the flow at the outlet of the brachial artery, R 2 is the flow resistance at the outlet of the brachial artery, L 2 is the inductance or inertia at the outlet of the brachial artery, v out is the pressure at the outlet of the brachial artery, C p is the peripheral elasticity, L p is the peripheral inductance, and v R is the peripheral pressure. These subscripts are also identified in Figures 1-3, depending upon its application. The model equations for fourth order WK model ( Figure 3) are: The model equation for the third order WK model ( Figure 4) is: For laminar flow, Poiseuille's law is used, given by where ΔP is the pressure drop along the artery, R is the resistance, is the fluid kinematic viscosity, is the fluid density, L A is the length of the artery, Q is the flow rate in the artery, and r is the artery radius. For turbulent flow, Blasius' equation 12 is used, given by, To include an entry flow effect Shah's equation for entry flow is used, given by where, the parameters h 1 through h 5 are defined as: In this equation, D A is the diameter of the artery, and c is a dimensionless constant equal to 0.000212. 4 When this model is used, the group (i 1 R 1 ) in Equation (1) is replaced with The group (i 2 R 2 ) in Equation (3) is altered for the different flow equations.
The capacitance values are obtained from the physiological volume change caused by a small pressure change: where Vol is the volume of blood in the given component. The brachial artery inductance values are obtained from An expression for L p is less straightforward because it involves multiple branching. For the sake of this study, a value of L p = L 1 was chosen. Simulations were also run with other values of L p , and even a change by 5 orders of magnitude, did not affect the model outputs. 4 The lumped-parameter equations were solved with the built in ODE45 and ODE113 functions in MATLAB. 24 Butcher's explicit sixth order RK method was manually coded in MATLAB. The models ran in approximately 30 s for each 5 s flow waveform.
Butcher's explicit sixth order RK method equations are given by: and where k1, k2, k3, k4, k5, k6, and k7 are the slopes, h is the step size, and f(x,y) is the function. The subroutine "findpeak" is built in MATLAB and used to obtain the number of oscillation peaks in the flow waveform. It returns a vector with local maxima (peaks). A local peak is a data sample that is either larger than its two neighboring samples or equal to infinity. 24 Error percentage is calculated as: where, Q ave is the average flow rate.

Part-I Comparison with CFD result
In the presented flow and pressure waveforms, the transient part (first flow cycle) is excluded.

Inlet pressure waveforms
The inlet pressure waveform is obtained as a reconstructed waveform from the first four harmonics from Stergiopulos et al. 4,13 Figure 5 shows the inlet pressure waveform that is analogous to the inlet brachial artery pressure with the systolic and diastolic pressure in the physiological range, that is, 125/78 mmHg.
F I G U R E 5 Inlet pressure wave (excluding transient)

3.1.2
Comparison between different orders of WK models using ODE45, ODE113, and Butchers' sixth order method In this section, brachial artery resistance is modeled by Poiseuille's flow. Figure 6 shows the flow rate at the outlet of the brachial artery for normal peripheral resistance. The flow rate generated by the fourth order WK model oscillates over time with a frequency of about 20 Hz ( Figure 6). The third order model has highest oscillations. The fifth order model exhibited an oscillation at about 5 Hz, but these oscillations correspond to similar features in the original pressure waveform. For low R p (Figure 7), the third order model again exhibits highest small oscillations, while the fourth order model exhibited low amplitude oscillations, with occasional spike-like anomalies. The fifth order model had the largest amount of oscillations, near 30 Hz. For both low and high R p , the flow rate for the fourth order model straddles that for the third order model. The amplitude of the fifth order model flow rate was larger than that of the third order model at low R p , but about the same at high R p . The maximum, minimum, and mean flow rates were smaller for the fifth order model than for the other two models at low R p .
Flow rate at the outlet end of the brachial artery for normal and low peripheral resistance using ODE113 is shown in Figures 8 and 9, respectively. For low R p , the third and fourth order models led to significant noise in the flow waveforms, and for high R p , all models led to flow waveform noise, with the fifth order model noise somewhat lower than that of the other two models. The third and fourth order models generally tracked one another. For normal R p , the fifth order model led the other two models in time and again exhibited the 5 Hz waves during deceleration and a lower minimum velocity. For high R p , the fifth order model generally tracked about 35ml/min lower than the other two. Peak and mean Reynolds numbers werearound 200 and 150 for low R p , and they were around 1500 and 1150 for low R p .
When Butchers' sixth order RK algorithm was used, the third and fourth order models were indistinguishable from one another (Figures 10 and 11). The noise and oscillations seen previously were not present, except for the 5 Hz waves during deceleration. These waves were most prominent for the fifth order model at normal R p and were about equal for all three models for low R p . The fifth order model led to a larger amplitude at normal R p while at low R p , amplitudes for all model orders were about the same, but the mean value was lower for the fifth order model.To reach conversion at lower R p , step size was reduced from the value of 0.001 s used in all other simulations to 0.0008 s.

Effect of step size on the ODE solutions
The effect of step size on the ODE solutions was examined for the fifth order WK model at low R p . The flow conditions were modeled by Poiseuilleflow. Brachial artery parameters are same as mentioned in Part-I. "NC" is No Convergence. Table 1 shows that ODE45 converges and gives result in most of the cases. Only h = 0.1 s led to a smooth curve. The other step sizes led to oscillations in the waveform at approximately 28 Hz with an amplitude of about 40 ml/min. Table 2 shows results from Butcher's sixth order method. The method converged for step sizes greater than 0.01 s, and the waveforms for all converged step sizes were nearly identical. The peaks were correlated with inflection points in the input pressure waveform and are thus considered to reflect the flow physiology. Table 3 shows results from the ODE113 solver. The method converged for step sizes greater than 0.1 s, but the waveforms for steps of 0.01 s and greater show fluctuations with a random character that are as large as 80 ml/min and uncorrelated with the pressure waveform. These results show that ODE113 is not a reliable option for this kind of problem. Step size (h) P max (ml/min) P min (ml/min) M (ml/min) Peaks in first cycle Peaks in three cycles

Wall displacement
In the previous study, 4 the radial wall displacement at the brachial artery inlet, obtained from the fifth order, WK model (with ODE45 as a solver) for normal R p at the inlet region closely matched that obtained from the CFD simulation. Therefore, this section compares the wall displacements obtained from different ODE solvers. Figure 12(A) shows radial wall displacement at the brachial artery inlet from the fifth order model for three ODE solvers. The ODE45, ODE113, and Butcher sixth order solvers exhibited almost no difference.

TA B L E 3
Features obtained from the ODE113 solver for different step sizes with low R p Step size (h) P max (ml/min) P min (ml/min) M (ml/min) Peaks in first cycle Peaks in three cycles at the brachial artery inlet from the ODE45 solver for third, fourth, and fifth order WK models for normal R p . (C) Wall displacement at the brachial artery inlet from the ODE113 solver for third, fourth, and fifth order WK models for normal R p . (D) Wall displacement at the brachial artery inlet from Butchers sixth order solver for third, fourth and fifth order WK models for normal R p Figure 12(B) shows the radial wall displacement of the brachial artery from the third, fourth, and fifth order models, solved with ODE45. The fourth order model produced strong oscillations at approximately 20 Hz. The third order model produced smaller oscillations at a much higher frequency that appear to occur in bursts. The fifth order model showed no discernable oscillations. Figure 12(C) shows the radial displacement at the inlet for different WK models using ODE113. The third order model exhibits the strongest fluctuations, and these have a random appearance. The fourth order model also led to fluctuations, but they were at a much lower frequency. The fifth order WK model shows no discernable oscillations. Figure 12(D) shows the radial displacement at the inlet for different WK models using Butchers method. All of the WK models generated similar results with less oscillation. However, the displacement predicted by fifth order model was slightly different, particularly in the diastolic region. Figure 13(A) shows radial displacement at the mid-artery in the fifth order WK model predicted by different solvers for low R p . In contrast to normal R p (Figure 12(A)), the ODE45 and ODE113 solvers show substantial oscillations. Figure 13(B) compares the radial displacements of the inlet and mid-artery for low R p in fifth order WK model using Butcher's sixth order method. Because pressure decreases with distance along the artery, the wall displacement is smaller at the mid-artery. This result is also relevant to the CFD results. 4 Figure 13(C) shows radial displacement from the fifth order WK model for normal R p , as solved with Butcher's sixth order method in inlet and mid artery. The inlet and mid artery displacements are nearly identical for normal R p because F I G U R E 13 (A) Radial displacement at the mid-artery in fifth order Windkessel model for different solvers (low R p ). (B) Radial displacement for low R p in fifth order Windkessel model using Butchers sixth order (Inlet vs Mid artery). (C) Radial displacement for normal R p in fifth order using Butchers sixth order (Inlet vs Mid-Artery) the peripheral resistance, rather than the brachial artery itself is the main source of pressure loss. This result agrees with the CFD results. 4

Part-II Comparison of numerical and clinical results
Q max , Q min , and Q ave are the maximum (systolic) minimum (diastolic), and average flow rates between 1.95 s and 5 s, respectively. The outlet of the brachial artery, which is connected with the anastomosis, can have a different flow pattern than the inlet. Therefore, the brachial artery was divided into an inlet region 13.3 cm long and an outlet region 6.67 cm long. Although irregular blood flow patterns are seen within the anastomosis the blood flow in the outlet region can be laminar or turbulent, and which regime exists can be difficult to predict at Reynolds numbers straddling 2000. 5

Comparison of the flow rate obtained with 0D lumped parameter model versus clinical data
Caroli et al reported data for 63 patients undergoing hemodialysis. 7 Clinical data for the brachio-cephalic side to end (BCSE) AVF, shown in table 2 from Caroli et al, 7 are used for the comparison. For inlet pressure, diastolic and systolic pressure from table 1 of the same paper is used. 7 Table 4 shows different brachial artery radii and the assumed peripheral resistances that were used to simulate the flow. The blood pressure waveforms used at the inlet are shown in Figure 14. Condition A models the brachial artery 1 day post operation, and the other models conditions 40 days post operation. 7 The brachial artery pressure for A is assumed to be approximately 94/60 mmHg (Systolic/Diastolic). The other blood pressure parameters are approximately equal to those reported in Caroli's paper. 7 The peripheral resistances given in Table 4 are the values required to obtain the reported flow rate. Table 5 reports the flow rates for conditions A through D with brachial artery input and output modeled as Poiseuille flow. For low blood pressure and smaller brachial artery radius, the values of flow rates generated by the simulation differ from clinical data by higher percentage. However, when inlet pressure and arterial radius are high, the percent difference is less. Figure 15 shows the flow waveforms. Table 6 reports the flow rates for conditions A through D with brachial artery inlet region modeled by Shah's entry flow equation and outlet region modeled by Poiseuille flow equation. Again, the pre-op conditions (D) led to a closer match between the measurements and the simulation. The difference in the flow rate for low pressure condition was slightly less than when both inlet and outlet were modeled with Poiseuille flow. Figure 16 shows the flow waveforms.    Table 7 reports the flow rates for conditions A through D, with brachial artery inlet region modeled by Shah's entry flow equation and outlet region modeled by Blasius flow. The error percent has reduced significantly for low blood pressure conditions. The error percent was higher than the previous values only for condition D. Figure 17 shows the flow waveforms. Table 8 reports the flow rates for conditions A through D with brachial artery inlet region modeled by Poiseuille flow and outlet region modeled by Blasius flow. The percent error has significantly decreased for conditions C and D. However, it has slightly increased for conditions A and B as compared with Table 7. Figure 18 shows the flow waveforms.

DISCUSSION
The fifth order WK model solved with Butcher's method to model the outlet flow rate of the brachial artery in an AVF generated least oscillation and most consistent flow rate of the MATLAB ODE solvers. While the MATLAB ODE solvers have better convergence range than the manual Butcher's method, can adjust the step size automatically, and are more user friendly, the fluctuations in the average flow rate and oscillation peaks (in the flow waveform) are drawbacks of ODE45 and ODE113 solvers. MATLAB ODE solvers also generated inconsistent flow rate with varying step sizes. Decreasing the step size 'h' did not generate consistent systolic and diastolic flow rate or reduce the oscillation peaks. Manual Butcher's method has a fixed convergence range, however once it converges it will generate consistent flow rates. The result from the Butchers' method has been validated with the CFD results 4 for normal and low R p . The total number of peaks between 2 and 5 s should be three times the number of peaks between 2 and 3 s if the flow solution is consistent from cycle to cycle. However, this ratio was not obtained, indicating that the MATLAB ODE solver is generating inconsistent changes in the oscillation pattern per flow waveform. Butcher's method, however, is consistent in the flow waveform pattern. The average flow in all the three solvers is similar. Manual Butchers' method requires trial and error to adjust the step sizes until it reaches convergence. ODE113 is not recommended for this kind of problem. Physiologically, spikes should not appear in the flow unless flow disturbances are present, as the shape of the flow waveform also signifies the various physiological conditions. 14,15 The source of the oscillations is still unknown and has to be resolved by the MATLAB software developers themselves.

4.1
Accuracy of the fifth order WK model

Normal peripheral resistance
Third and fourth order WK models solved with MATLAB ODE solvers for normal R p generated more oscillations than fifth order model. Peak flow rate predicted by the fifth order WK model is less than that of the fourth order model and higher than third order model. Oscillations generated by the third order WK model are higher than fourth and fifth order models. However, the average flow rates predicted by all of the three models are similar. Oscillations generated by ODE 113 for all the three models (using all the three solvers) are highest. Peak flow rate is also higher for third order WK model using ODE 113. It is interesting to note that the third and fifth order models showed closer values for ODE45 and Butcher's sixth order method. The fifth order model generated consistent peak and average flow rates, along with stable radial wall displacements, for all three solvers.

Low peripheral resistance
The oscillation peaks are larger in fifth order WK model for ODE113 and ODE45.The average flow rates predicted by third, fourth, and fifth order WK models solved withODE45, ODE113, and Butchers' method are 899, 901, and 864 ml/min, respectively. The fourth order WK model solved with ODE113 showed lower frequency, higher amplitude peaks, and higher peak flow Reynolds number. Generally, higher Reynolds number promotes turbulent flow, causing higher wall shear stress and lower blood flow velocity for a given flow rate. However, a definitive Reynolds number for transition cannot be established, as shown from this study, and caution has to be taken while selecting the solver and the WK model. In the previous study, 4 Shah flow in the lumped parameter model matched the solution from CFD in terms of peak flow rate than did Poiseuille flow for low R p . Shah has lower peak flow rate than Poiseuille when solved with ODE 45. In this study, with the Poiseuille model, the peak flow rate was 1085 ml/min from the ODE45 solution and 1059 ml/min from the Butcher's solution, demonstrating the influence of ODE solver on the accuracy of the flow rate.
Flow resistance influences the overall flow rate more strongly than inductance and capacitance and depends on the radius to the fourth power. Peripheral Inductance (L p ) and peripheral capacitance (C p ) havesmall influence on the flow rate, but the effect of capacitance on the flow waveforms can be seen for C p less than 10 −11 m 3 Pa −1 .

Impact of flow equations in the arterial segment
Caroli et al. suggested that flow rate less than 400 ml/min indicated a failed fistula in four patients. 7 Shah inlet with Blasius outlet conditions generated smallest outlet flow rate in the brachial artery for low inlet pressure, which suggests that turbulence at the end of the brachial artery proximal to the stenosis reduces blood flow in the anastomosis. Solution to conditions one day post-operation 7 (condition A) shows that low blood pressure and high peripheral resistance right after the creation of an AVF, with the small arterial diameter, decreases outlet flow rate. Localized blood pressure at the end of the brachial artery might also be an important factor to explain the performance of an AVF. With larger radius and higher blood pressure, Poiseuille flow conditions seem to predict the flow value more correctly. Vessel wall elasticity is important for the remodeling of the arteries. If the blood pressure is lower than normal, shear stress on the brachial artery wall will be smaller, which might impair the arteries capacity to dilate and ultimately might cause fistula failure. This study further examines the claim that entry flow modeled with Shah's equation is clinically relevant to the steal syndrome. 4 The solution of the Shah inlet and Blasius outlet along the varying length of the brachial artery plausibly models entry flow-related shear stress and wall shear stress with low R p is correctly modeled by Shah equation. 4 The high shear stress then leads to a low pressure, and therefore low flow rate, at the arteries feeding the hand.
This model suggests that the Poiseuille inlet and outlet and the Shah inlet with Poiseuille outlet (laminar flow equations) are not suitable for certain inlet pressure conditions or for smaller brachial artery radius. Outlet flow modeled with Blasius equations more closely matches the clinical findings. Therefore, this model showed that the inlet region of the brachial artery blood flow can be modeled with laminar flow equations while the outlet region can be modeled with turbulent flow equation. When the brachial artery was equally divided into inlet and outlet regions (result not shown) the flow rate did not changed for the Poiseuille inlet and outlet and the Shah inlet with Poiseuille outlet, but it changed significantly for Shah inlet with Blasius outlet and Poiseuille inlet with Blasius outlet. That is, turbulence proximal to the stenosis noticeably affected the flow rate. Arterial elasticity (capacitance) has small effect at higher flows.

Limitations on assumptions made in Part-II
Peripheral resistance is used as an estimated value to obtain AVF flow. 4 The 0D lumped parameter model is said to give accurate results only for diastolic pressure, where the wave propagation phenomenon is ignored, not necessarily for systolic pressure. However, the 0D model is still useful to model boundary conditions in terminal arteries, which represent blood microcirculation. The systolic and diastolic range is an assumed value based on the clinical data rather than an exact value.

4.4
Comparison with other studies

MATLAB ODE solvers
Abelman et al noted that ODE45 and other solvers in MATLAB are limited for large step sizes. 25 Table 1 shows the performance of the ODE solver in various step-sizes. In this study, ODE45 generated varying flows upon changing the step-size. Dallas et al compared the MATLAB ODE solvers for hybrid problems, that is, stiff and non-stiff problems. 18 Their findings differ with the findings of this paper because this paper deals with non-stiff problems. In contrast to their solver rankings, ODE113 performed poorly in the present study.
Fazio described the instant failure of ODE45 in solving a system of equations that consisted of two coupled oscillators. 26 In our study, ODE45 does not fail to provide a solution but did not give a correct solution and might provide misleading information.
Kulikova et al stated that the MATLAB ODE solvers do not take into account special features of differential equations such as symmetry and semi positivity of the error covariance matrix. 27 They further mentioned that specially developed solvers outperform MATLAB ODE solvers.

WK model order
Stergiopulos et al 13 analyzed the efficiency of third and fourth order WK models. Their study is more realistic as they have compared the results from the model with in vivo canine data and human aortic pressure waves. They concluded that the fourth order WK model matches the in vivo data more closely than the third order model. They used transmission line models and the concept of pulse wave propagation theory, which was incorporated in the third and fourth element WK model. Hauser et al 28 combined different elements of the WK model with MATLAB Simulink to study blood flow. The authors concluded that the most realistic blood pressure curves were obtained from the fourth order WK model. Their model was based on flow rate as an input and pressure as an output. They explicitly state the resistance values, the model's physiological application, or how the model was derived. However, this study focuses on the AVF and the blood flow resistance. Other flow parameters are adjusted to match those of an AVF. Kokalari et al 11 described various WK models used in diagnosis of cardiovascular diseases. They noted the efficiency of the fourth order WK model over third order WK model. Furthermore, they indicated a 5-element WK model as better describing the dynamics of blood microcirculation.

CONCLUSION
Third, fourth and fifth order lumped parameter models were used to simulate the brachial artery blood flow in an AVF. ODEs of the model were solved using ODE45, ODE113 and Butcher's sixth order RK method in MATLAB.
Poiseuille's law, Shah's law, and Blasius's equation were used to define the blood flow resistance in the brachial artery. Shah inlet with Blasius output flow conditions along the varying length of the brachial artery generated smaller flow rates, which may explains maturation failure of an AVF for the case presented in this paper. Butcher's sixth order method showed stable, consistent and optimal oscillation in the flow waveform as compared with the MATLAB built in solvers, ODE45 and ODE113.One reason for the variations is that a higher order RK method employs slopes at more end points of the interval of the extrapolation equation. Despite the stable results, Butcher's sixth order RK method requires trial and error to reach convergence. Physiological oscillations or vibrations in an AVF flow are subject of future research. However, they should not arise from a limitation of the ODE solver but from some anatomical abnormalities in an AVF itself. Third and fourth order lumped parameter models generate similar result only at normal peripheral resistance when compared with a fifth order model.
In this paper, average blood flow rate was obtained with various blood flow resistance models and compared with clinical results. It was found that a segment of the brachial artery can be modeled with 0D lumped parameter model, with some adjustments. This study showed that flow resistance equations also play significant roles in modeling the blood flow in brachial artery. This model showed that brachial artery dilatation is crucial in determining the successful fistula. According to the findings of this simulation smaller brachial artery radius, higher peripheral resistance, low blood pressure and turbulent flow at the outlet end of the brachial artery near the anastomosis are indicators of a failing to mature AVF. Therefore, it is subject of further research to understand the relationship between arterial dilatation, blood pressure, and peripheral resistance.

PEER REVIEW INFORMATION
Engineering Reports thanks Raheem Gul and other anonymous reviewers for their contribution to the peer review of this work.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CONFLICT OF INTEREST
The authors have no conflict of interest relevant to this article.