- 1Key Laboratory of the Ministry of Education for Optoelectronic Measurement Technology and Instrument, Beijing Information Science and Technology University, Beijing, China
- 2Guangxi Key Laboratory of Regenerative Medicine, Research Centre for Regenerative Medicine, Guangxi Medical University, Guangxi, China
- 3School of Mathematics and Statistics, University of Glasgow, Glasgow, United Kingdom
- 4Second Affiliated Hospital of Nanjing Medical University, Nanjing, China
One-dimensional (1D) hemodynamic models of arteries have increasingly been applied to coronary circulation. In this study, we have adopted flow and pressure profiles in Olufsen's 1D structured tree as coronary boundary conditions, with terminals coupled to the dynamic pressure feedback resulting from the intra-myocardial stress because of ventricular contraction. We model a trifurcation structure of the example coronary tree as two adjacent bifurcations. The estimated results of blood pressure and flow rate from our simulation agree well with the clinical measurements and published data. Furthermore, the 1D model enables us to use wave intensity analysis to simulate blood flow in the developed coronary model. Six characteristic waves are observed in both left and right coronary flows, though the waves' magnitudes differ from each other. We study the effects of arterial wall stiffness on coronary blood flow in the left circumflex artery (LCX). Different diseased cases indicate that distinct pathological reactions of the cardiovascular system can be better distinguished through Wave Intensity analysis, which shows agreement with clinical observations. Finally, the feedback pressure in terminal vessels and measurement deviation are also investigated by changing parameters in the LCX. We find that larger feedback pressure increases the backward wave and decreases the forward one. Although simplified, this 1D model provides new insight into coronary hemodynamics in healthy and diseased conditions. We believe that this approach offers reference resources for studies on coronary circulation disease diagnosis, treatment and simulation.
1. Introduction
With the ever increasing load of coronary heart disease, studying the coronary arteries is of importance. Coronary blood flow and pressure are key components in cardiac pump function. Clinical measurements, particularly non-invasive ones, are difficult to perform and necessarily limited in terms of resolution (Toyota et al., 2005; Chiribiri et al., 2013). Invasive dynamic measurements are regularly performed in the coronary circulation system, and in many cases may be clinically indicated, but they do carry some risk, so non-invasive approaches would be preferable in circumstances where invasive measures are not indicated. Mathematical modeling of coronary vessels is an alternative approach to understanding the relationship between coronary flow/pressure and physical properties of the cardiovascular system, at required locations, and it also provides useful insights, such as feedback control of the coronary vascular resistance and wave intensity.
Over the last twenty years, various mathematical models have been developed for cardiovascular problems, ranging from one segment of the vessel to the whole vessel system (Pietrabissa et al., 1996; Kim et al., 2010; Mynard, 2011; Karimi et al., 2013). Although most of the models are developed for systemic circulations, the mathematical principles are similar for coronary circulations. The simplest approach is the so called lumped-parameter or zero-dimensional (0D) modeling. This is based on a resistance-compliance Windkessel approach for the circulation system (Frank, 1899; Pietrabissa et al., 1996). This type of model has been used to simulate the cardiovascular system and the distribution of the blood flow and as an investigative tool to evaluate the different bypass treatments (Pietrabissa et al., 1996). Zero-dimensional models are also widely adopted as out-flow boundary conditions to provide a physiologically accurate pressure or flow profile (Huo and Kassab, 2007; Alastruey et al., 2008; Mynard, 2011). It is particularly useful for evaluating the effects of micro-circulation in patient-specific simulations (Alastruey et al., 2008). Many studies have also demonstrated that the myocardial contraction provides extra forces to the terminal coronary wall (Spaan et al., 2000; Algranati et al., 2010). Mynard (2011) established a comprehensive hemodynamic model for coronary arteries in which a 0D model was used at the boundary to reflect the effects of micro-vascular beds. Duanmu et al. (2018) recently developed a patient-specific 0D model of coronary circulation in which the root impedance was evaluated based on a static structured tree model, and good agreement with experiments was obtained.
The major limitation of the 0D approach is over-simplification in omitting the wave information. In addition, very few studies have considered the interaction between the ventricular contraction and the coronary arteries when modeling the coronary flow (Van der Horst et al., 2013; Mynard et al., 2014; Mynard and Smolich, 2015). This is because 0D models are unable to simulate the wave effects due to heart fluctuations (Van der Horst et al., 2013). Most of the coronary models use 0D terminal boundary conditions as feedback from ventricular contraction (Mynard et al., 2014; Mynard and Smolich, 2015), which is a crude approximation for the physiological vascular bed.
In one-dimensional (1D) models, the coronary arteries are treated as tapered tubes. The governing equations for the flow are derived from the axisymmetric Navier–Stokes equations, and the vessel walls are assumed to be either rigid or compliant (Olufsen, 1999). Owing to its simplicity and a yet more rational description of the physiology, 1D modeling has been widely used in circulation systems, such as the systemic arteries (Olufsen, 1999; Olufsen et al., 2000), microcirculation (Alastruey et al., 2008) and coupled models Huo and Kassab (2007); Mynard (2011). Olufsen et al. studied such a model for systemic circulation in which a structured tree boundary condition was used (Olufsen, 1999; Olufsen et al., 2000). A 1D model was also developed for subject-specific human pulmonary circulation (Qureshi and Hill, 2015). Similar approaches were used to predict flow and pressure in the pulmonary arteries (Olufsen et al., 2012), which showed that a decrease of the compliance in the large arteries leads to a pressure rise. Stergiopulos et al. (1992) showed that a 1D model can provide a better result than the 0D model in accounting for the nonlinear properties of the arterial wall. 1D modeling is often used to build a coupled circulation system with frequency outlet condition. Recently, Chen et al. (2016) built a coupled system between a three-dimensional finite-strain left ventricle and a physiologically-based 1D model of systemic circulation, and studied the effects of rarefaction and arterial stiffness on the performance of the heart.
With the development of medical imaging technologies, detailed 3D coronary arterial models can now be reconstructed from in vivo imaging (Mittal et al., 2005; Raff et al., 2009). Simulations of 3D circulation systems have been performed with different outlet boundary conditions (Vignon-Clementel et al., 2006; Kim et al., 2010). Such an approach, however, is extremely time consuming and the accuracy of the result depends on the parametric boundary conditions, which are usually difficult to obtain. For this reason it is very important to continue to develop more accurate patient-specific 1D models for the coronary system, in order to provide a real-time clinical support system.
One of the advantages of 1D modeling is to be able to analyze wave intensity (WI), which is important in evaluating the myocardial and hemodynamic functions. The waveform of the pressure or flow in the arteries can provide a convenient way to assess the wave energy transmission in a time domain-based manner in the cardiovascular system (Parker, 2009). In addition, biventricular pacing that increases the coronary blood flow velocity from diastolic-dominant backward decompression waves can be detected though WI analysis (Kyriacou et al., 2012). Combined analysis of pressure and flow velocity using WI analysis demonstrates that coronary hemodynamics are naturally most reflective of the upstream fluid dynamics of a stenosis Sen (2013). However, WI analysis based on measured flow and pressure alone can sometimes be inaccurate or even misleading (Sun et al., 2000; Sen et al., 2014; Mynard and Smolich, 2015) due to the difficulties in measuring coronary flow and pressure in vivo. When used in a robust computational model, WIA can quantify contributions of pressure and flow waveforms at arbitrary locations and different physiological conditions (Willemet and Alastruey, 2015).
In this paper, we have developed a 1D model of the large coronary arteries based on in vivo CT scans. This model takes account of the detailed coronary geometry and the pulsatile pressure loading within the myocardium. Unlike previous coronary models with boundary conditions estimated from 0D models (Mynard et al., 2014; Mynard and Smolich, 2015), we use a structured-tree model (Olufsen, 1999; Olufsen et al., 2000) for the vascular beds of smaller coronary vessels to provide the boundary conditions relating pressure and flow at the distal end of each terminal artery, added to an additional pressure feedback term resulting from the ventricular wall contraction. We perform WI analysis of the left circumflex artery under different pathological conditions, such as increased arterial wall stiffness and vascular rarefaction in the structured trees. Finally, the effects of the feedback pressure and geometry uncertainties are investigated.
2. Methods
2.1. Reconstruction of the Coronary Vessels
The right and left coronary trees are reconstructed from the in vivo data from a male patient (age 61, 80 kg). GE 16-slice CT angiography was performed and medical imaging software (GE, AW) was used for visualization and measuring vessel length and radius. Blood flow data for the left circumflex branch of the coronary vessels was measured by an intra-coronary Doppler guide wire (FloWire, 1400-Floppy). The measurements were performed by The Second Affiliated Hospital of Nanjing Medical University, China. All the procedures were conducted according to established ethical guidelines (Raff et al., 2009). This study was approved by the Human Experiment Committee of The Second Affiliated Hospital of Nanjing Medical University. The patient provided written consent.
In order to reconstruct the coronary tree structure, each vessel was rebuilt along the center line from the CT angiography data at diastole during maximum hyperemia. The coronary tree is separated into segments according to the clinical anatomy following the American Heart Association Standard (Raff et al., 2009). The coronary trees are separated into two branches: the left main coronary artery (LMCA), which includes the left anterior descending branch (LAD), the left circumflex branch (LCX), the diagonal artery (DIAG), and the marginal branch (MARG), and the right coronary artery (RCA) which includes the posterior descending artery (PDA) and the posterior lateral artery (PLA). In the left branch, each segment is further divided into several segments as shown in Figure 1. Each vessel is tapered from the proximal to distal ends. The measured geometries of all the vessels are listed in Table 1.
Figure 1. The structure of the coronary tree, including the right and left coronary arteries Raff et al. (2009): (A) CT reconstruction of coronary circulation, (B) the 1D coronary tree model, (C) modeling of the trifurcation, and (D) the terminal branches are connected to a capillary bed in the myocardium using a structured tree model with feedback pressure pf(t).
Following Olufsen (1999) and Olufsen et al. (2000), the distal end of each terminal large artery is connected to a vascular bed (represented in the model by a structured tree) of smaller vessels that lie within the myocardium. In the structured tree, each parent vessel of radius rp is divided into two daughter arteries with radii rd1 = αrp and rd2 = βrp, where 0 < β < α <1. The branching process continues until the radii of daughter vessels reach the minimum radius, rm. The bifurcation process obeys:
where ξ is the radius exponent, η is the area ratio, and γ is the asymmetry ratio. It should be noted that the parameters are related.
2.2. 1D Coronary Model
In the 1D model, we consider the vessel wall to be a deformable tapering tube, whose radius and wall properties are functions of the centerline distance x from the proximal end of the vessel. We further assume the vessel wall is impermeable, and that blood is an incompressible Newtonian fluid. The continuity and momentum equations are:
where A is the cross-sectional area, q = uA is the flow rate, u is the velocity, ρ is the blood density, μ is the blood viscosity, r(x, t) is the vessel radius, and δ is the thickness of the boundary layer (Olufsen, 1999; Olufsen et al., 2000). The blood pressure p is related to A in the tube-law:
where r0 is the vessel radius, h is the wall thickness, p0 is the reference pressure, and A0 is the cross sectional area when p = p0, and E is the Young's modulus, which is a function of r0 given by:
where the coefficients k1, k2, k3 are taken from Olufsen (1999). These equations are solved numerically using a Lax–Wendroff method (Lax and Wendroff, 1960).
Unlike the structured tree models for systemic and pulmonary circulations, where only bifurcations exist, a trifurcation is observed in the coronary circulation at the end of LMCA (consisting of LCX, LAD and DIAG) as shown in Figure 1C. To use the same branching process (Equations 1, 2), we divide the trifurcation into two bifurcations. Specifically, the LMCA is separated into daughter 1 (the DIAG) and daughter 2 with a radius of 0.185 cm derived from (1). The daughter 2 is further divided into two granddaughter vessels: LCA and LCX. The length of daughter 2 is small (0.25 cm) so that the structure is close to the original trifurcation.
2.3. Impedance of the Vascular Beds
The structured-tree model is used as the outflow boundary of each terminal large vessel as shown in Figure 1D. In the frequency domain, the linearized momentum and continuity equations for the small arteries in the terminal vascular bed are (Olufsen, 1999):
where ω is the frequency and K is a parameter related to the Womersley number (Olufsen, 1999), is the vessel compliance, and P(x, ω), Q(x, ω) are the complex Fourier transforms of p(x, t) and q(x, t), and ω takes the values that are integer multiples of the cardiac frequency.
From (7) and (8) we obtain:
where c is the wave propagation velocity in the small vessels of the structured-tree,
The impedance Z(x, ω) is computed from:
At a vessel bifurcation, conservation of flow and continuity in pressure lead to:
where Zp is the impedance of the parent vessel, and Zd1 and Zd2 are the impedances of the daughter vessels of diameters d1 and d2. Using (9) and (12), we can deduce the root impedance for a vessel of length l:
where Z(l, ω) is the impedance at the distal end of the vessel, and . Thus, (9), (12) and (13) can be solved recursively to obtain the root impedance of the terminal vessel of the vascular beds (Olufsen et al., 2000).
2.4. Inlet and Boundary Conditions
For a given ventricular pressure, the inlet flow rates to the left and right coronary branches are estimated from a lumped parameter model, for a heart rate of tc = 0.9s (Pietrabissa et al., 1996). The coronary blood flow is 4–5% of the cardiac output (Hall, 2015). We determine the left ventricular and right ventricular pressures, and following (Avanzolini et al., 1988):
where is the peak isovolumetric pressure at the volume , is the time-varying elastance, Ed is the diastolic passive elastance, is the systolic passive elastance, the activation function σ(t) = [1−cos(2πt/ts)] during systole and 0 during diastole, ts = 0.4s is the time at end of systole, V is the ventricular volume, is the rate of volume change, and is the resistance of the ventricular myocardium. A similar expression holds for .
The reference pressure p0 in (5) is set to the time-averaged value of the left or right ventricular pressure for the left and right coronary arteries, respectively (Mynard et al., 2014). The large coronary arteries in this model mostly lie outside the myocardium. The terminal large arteries each connect to a vascular bed of smaller vessels represented by a structured tree. The contraction of the myocardium during systole squeezes the vascular bed and thus increases the pressure at the end of each terminal artery, where it meets its vascular bed. To account for the effects of the myocardial contraction on the coronary vessels, we introduce a time-dependent feedback pressure pf(x, t) that is added to the distal boundary conditions for each terminal large artery, i.e., where each terminal vessel is connected to its vascular bed (represented by a structured tree); this increases the resistance to flow due to the compression of the vascular beds during systole (Mynard, 2011). We choose pf to be proportional to ventricular pressure plus the coronary bed pressure ps (Mynard and Smolich, 2016), i.e.,
where ps = 20mmHg (Hoffman and Spaan, 1990). The ventricular-pressure feedback ratio φ depends on the location of the vascular bed supplied by the terminal vessel because the left and right ventricle peak-systolic pressures are different; see Table 2. The chosen values give physiologically realistic pressure profiles.
Table 2. The ventricular-pressure feedback ratio is applied to all terminal coronary vessels, and depends on the location of the vascular beds within the myocardium. The values are chosen so that the resulting pressure profile is physiological.
2.5. The WI Analysis
The net WI is a synthesis of different waves, for example the WI may consist of a large forward wave minus a smaller backward one. The primary forward wave is caused by ventricular contraction and transmitted via the aorta into the coronary tree. Backward propagating waves are caused by reflections at bifurcations and also by myocardial contraction in the vascular beds.
To perform the WI analysis, the forward and backward wave fronts of the pressure and velocity are defined as in Qureshi and Hill (2015):
over the sample times ti = iδt ∀i ∈ [0, N], ti ∈ [0, T], where is the computed pulse wave velocity in the artery. The wave intensity is then (Qureshi and Hill, 2015)
Waves are assumed to be linearly composed, and the net effect of the waves can be studied by analyzing the forward compression and decompression waves (FCW, FDW), and the backward compression and decompression waves (BCW, BDW).
The physiological events that determine the WI patterns of the coronary arteries can be distinguished by six phases (Davies et al., 2006): (1) During early systole, the coronary WI pattern begins with a manifest BCW and decreasing LCX velocity. This BCW is generated by myocardial contraction, which compresses the small vessel beds (Chilian and Marcus, 1982). (2) The systolic FCW caused by contraction of the left ventricle moves from the aorta to the coronary vessels. (3) Later in systole, a secondary BCW is probably generated by reflection sites such as distal coronary artery pulsation or vessel branches. (4) In proto-diastole, the myocardium begins to relax and generates an FDW which lowers the pressures in the coronary arteries. (5) Concomitantly, a dominant BDW is thought to result from ventricular relaxation; the small vessels are decompressed, which results in lower pressure at the vascular bed and a lower resistance of the micro-circulation. (6) At the end of diastole, it is speculated that a late FCW is generated from the proximal end of the coronary arteries due to the closure of the heart valve. In identifying these wave components, WI has proved to be an efficient diagnostic tool.
3. Results
3.1. The Parameters
All the parameters and their values used in the model are listed in Table 3 along with the corresponding references.
3.2. The Baseline Coronary Flow and Pressure
Figures 2A,B show the computed profiles of the pressure and flow rate in different coronary vessels throughout the cardiac cycle. At the beginning of systole, despite the rapid increase in the pressure, there is a sharp dip in flow-rate curve of most of the vessels; this is caused by the active heart contraction. The flow rate then increases with the pressure until it peaks, then both drop at the end of systole at about 0.4 s. In diastole, although pressure remains low, the flow rate increases in most of the coronary vessels due to the myocardial relaxation and pf(t) drops. However, the flow rate in the smaller vessels of the coronary trees, such as LDA3, remains low at 0.5 ml s−1, presumably because they are less affected by the myocardial relaxation. In the right coronary vessels (e.g., RCA) the flow rate has a different trend as the pressures and contraction are lower in the right ventricle than in the left ventricle.
3.3. Comparison Between Experiments and Published Results
The simulated pressure and flow waveforms in the main vessels of coronary tree are compared with our own in vivo measurements in the LCX, as well as in vivo experiments (Duncker and Merkus, 2004; Spaan, 2012; Pijls and De Bruyne, 2013) and 0D (Duanmu et al., 2018), and a 3D Kim et al. (2010) waveform from healthy adults (Figure 3). The in vivo flows in LAD and RCA are from Figure 3 of Duncker and Merkus (2004), and the pressures in LAD and RCA are separately from Figure 5.13 and Figure 6.10 in Pijls and De Bruyne (2013). The cardiac cycles of the clinical data are adjusted to have the same heart beat of 0.9 s. The simulated flow rate agrees well with our own experimental data in the LCX. Although the geometries of the coronary structure are different, our model also compares reasonably well with the published experimental measurements.
Figure 3. Simulation results (black solid) for selected coronary vessels [The pressure in LAD (A), LCX (C), and RCA (E). The flow rate in LAD (B), LCX (D), and RCA (F)] are compared with published in vivo data (red solid dots), the corresponding 0D model (dashed dot blue), and a 3D model (dashed pink).
The flow profiles from our 1D simulated model are different from those of the 0D (Duanmu et al., 2018) and 3D (Kim et al., 2010) models in some vessels. As shown in Figure 3D, the 0D flow rate of LAD is different during systole, which might be because of the lumped parameter assumption. In right coronary vessel RCA, the flow rate from the 0D and 3D models becomes systolic dominant, which might result from the 0D outlet boundary condition, and is not consistent with the experimental data in Figure 3F. Overall, our model compares reasonably well with the published experimental measurements (Duncker and Merkus, 2004; Spaan, 2012; Pijls and De Bruyne, 2013).
Diastolic dominant coronary arterial flow is also captured in our model. It is caused by decompression of the terminal bed consequent to myocardial relaxation (Van der Horst et al., 2013), which compresses the terminal vascular bed. This myocardial contraction's effect on the coronary flow is incorporated in the model using the feedback pressure. The flow in the right coronary vessel is less dominant in diastole, which is due to the lower blood supply to the ventricle from RCA.
3.4. Wavefronts and WI Analysis
We now apply the wave intensity analysis to identify the forward and backward traveling waves. Since most of the published in vivo measurements and modeling results are for the LCX vessel, we use the results for LCX as a benchmark case. In Figures 4A,B, forward and backward velocity and pressure components are calculated from the pressure and velocity waveforms. The WI results for the LCX are shown in Figure 4C, where positive values indicate compression and the negative ones decompression; is the forward flow and is backward. Six major WI components can be identified in Figure 4D, which are consistent with previous studies (Davies et al., 2006). Figure 4E shows that our results are comparable with the measurements by Davies et al. (2006).
Figure 4. (A) Simulated velocity (U0) and (B) pressure (P0) waveforms (solid black), decomposed into forward (“+” dashed pink), and backward (“-” dashed-dot blue) running components at the midpoint in the LCX. The separated wave fronts and intensities are shown in (C) and (D) (shaded red), and (shaded blue). (E) The total simulated wave intensity (black line) is also compared with the experimental data (red line with dots) of Davies et al. (2006).
Next, we show our results at three different locations of the coronary tree (RCA, LCD, LMCA) in Figure 5. The blood pressure in the RCA (Figure 5D) is almost the same as the pressure in the LCX (Figure 4), but the flow in the RCA (max 22.5 cm s−1) is lower. The six different waves in the RCA are shown in Figure 5J, the WI are much less than those in the LCX, which is also reflected by the less flow fluctuation in the RCA. In the LMCA, the flow velocity increases with a maximum value of 36 cm s−1 (Figure 5B), while the pressure is similar to those in the LAD and LCX, as shown in Figure 5E. The increased flow velocity in the LMCA also corresponds to the increased WI compared to the values in the LCX, shown in Figure 5K. The flow, pressure and WI in the LAD are very close to the LCX flow and pressure, as shown in Figure 5.
Figure 5. Simulated waveforms for different coronary vessels; the right coronary vessel RCA (1st column): velocity (A), pressure (D), wave front (G), and intensities (J). LMCA (2nd column): velocity (B), pressure (E), wave front (H), and intensities (K). LAD (3rd column): velocity (C), pressure (F), wave front (I), and intensities (L). Simulated velocity (U0) and pressure (P0) waveforms (solid black), decomposed into forward (“+” dashed pink) and backward (“-” dashed-dot blue) running components at the midpoint in each vessel. The separated wave fronts and intensities, and (shaded red), and (shaded blue), are displayed at the bottom.
The area enclosed by each of the major waves in one heart beat represents the flow energy. The wave energy analysis helps to estimate the physiological reaction of heart failure if coupled with other information such as coronary flow or pressure (Kyriacou et al., 2012; Lee et al., 2016). The flow energy from our model is compared with experimental studies from Davies et al. (2006) and Lee et al. (2016) in Table 4. The waves (a), (b) and (c) in Table 4 are due to systole, and there is considerable variation between all sources, probably due to differences in the measurement location and vasculature structure (Lee et al., 2016). The dominant wave BDW (e) contains the greatest energy ration, which agrees well with the two experimental studies. The FDW is due to reflection at junctions between vessels.
3.5. Pathological Cases
We now compare the baseline results with the pathological cases. Four different pathological cases are considered. (A) Stiffer large arteries which are associated with a higher risk of cardiovascular diseases such as stroke and coronary atherosclerosis (Weber et al., 2004). This is modeled by increasing k3 by 100% in Equation (6). (B) Rarefaction (reduced vessel density) in the terminal vascular bed is modeled by reducing the exponent ξ from 2.76 to 2.4 in (1), corresponding to a 20% reduction in small vessels (Olufsen et al., 2012). In borderline hypertension and normal young offspring of parents with hypertension, it has been suggested that micro-circulatory rarefaction may be a cause of hypertension (Olufsen et al., 2012). (C) Changed feedback pressure is associated with heart pump function and can be indicative of several coronary epicardial diseases (Siebes et al., 2009). This is considered by varying φ from 0 to 60% for LAD1; in particular, when φ=22%, LAD1 is connected to septum (as shown by our CT scan), and when φ=33%, LAD1 is connected to the left ventricle which sometimes occurs (Mynard and Smolich, 2016). (D) Variation in coronary geometry can have an impact on the model results; the changes could be pathological or due to inaccurate measurements (Uus et al., 2015). This is assessed by shrinking or enlarging the arterial length and radius by ± 5 and ±10%, respectively.
The comparison with cases A and B for the LCX are shown in Figure 6. We can see that a higher stiffness in Case A results in an increased velocity (30–35 cm s−1 peak), but pressure changes are moderate. In Case B, pressure is significantly increased (from the peak value of 120 to 150 mmHg), but velocity changes are not significant. It is interesting to see more information is provided by WI analysis. In particular, the difference between the Baseline and Case B is smaller in both pressure and velocity waves, but more changes are shown in Case A.
Figure 6. Comparison of results along LCX are shown for the Baseline (left): velocity (A), pressure (D), wave front (G), and intensities (J). Case A (middle): velocity (B), pressure (E), wave front (H), and intensities (K). Case B (right): velocity (C), pressure (F), wave front (I), and intensities (L). Simulated velocity (U0) and pressure (P0) waveforms (solid black), decomposed into forward (“+” dashed pink) and backward (“-” dashed-dot blue) running components at the midpoint in the LCX of different conditions (1st and 2nd row). The separated wave fronts and intensities, and (shaded red), and (shaded blue), are displayed in the 3rd and 4th rows. The black line in the bottom panel is the total sum of the WI.
Figures 7A,B show the comparison of the results for LCX between the Baseline (φ = 0.33) and Case C (φ = 0.0−0.6). Clearly, the value of φ has a great impact on the simulated results. Figure 7 shows that the coronary flow increases as φ decreases during systole, and vice versa during diastole. The dynamic feedback pressure is systolic dominant, and the blood flow increases when the feedback pressure is lower.
Figure 7. Flow (A) and pressure (B) in the LCX for different values of φ: 33% when the LCX is connected to the ventricular and 22% when connected to the septum.
Finally, we compare the baseline with Case D, where the overall size of the coronary tree (including all the radii and lengths) is changed to (90, 95, 105, and 110%) of the baseline case. We observed that although the wave profiles are similar in all cases (not shown), the magnitudes are changed, as shown in Table 5. We first note that flow resistance is approximately proportional to vessel length. We now explain the influence of size on WI. The arterial stiffness (elastic modulus) is the inverse of compliance, and a stiffer vessel has a higher FDW (Davies et al., 2006). The decrease in size is equivalent to a compliance reduction which results in augmentation of stiffness. Hence, the increase/decrease of the geometrical measurement corresponds to a low/high FDW.
Table 5. Effects of LCX size changes. i.e., changes in length and distal radius. The maximum flow, pressure and WI are given as the percentages of the corresponding baseline values.
4. Discussion
In this study, the systemic arteries model originally developed by Olufsen (1999) has been extended to simulate the dynamics of the anatomically based coronary arteries. Although a number of other multiscale simulations of coronary flow have been carried out (Mantero et al., 1992; Mynard and Smolich, 2015), none of these have considered a structured-tree model as in Olufsen (1999) and incorporate myocardial contraction as feedback pressure in the vascular network. The pattern of WI in coronary arteries is quite different from the systemic circulation due to myocardial contraction, and the net WI of the LCX agrees well with experimental data, as shown in Figure 4E.
The simulated waveforms demonstrate that the flow velocity and WI patterns differ between the left and right coronary vessels. The shapes of the flow, pressure and wave intensity curves in the left coronary vessels (e.g., columns 2 and 3 in Figure 5) are the similar, but with different amplitudes. The greatest flow occurs during diastole. Larger flow and WI in the LMCA can be observed during diastole since it is close to the aortic root. Compared with left coronary vessels, the WI in the RCA is different and much smaller. This is because there are fewer structural branches and lower terminal pressure differences (Hadjiloizou et al., 2008) in the right coronaries.
We also model some pathological cases and show that a higher flow is always associated with an increase in the three forward wavefronts in systole. Our results also confirm that increased stiffness of the coronary arteries will increase the FDW during diastole, consistent with observations by Davies et al. (2006). The stiffened arteries can create an energy surge into coronary vessels. Hence, the FDW might be a good indicator for diagnosis of diseased coronary arteries.
The augmentation of pressure can be seen in rarefaction even though there is no obvious change in the flow. However, rarefaction contributes to hypertension with the increased BCW and BDW. This abnormally high pressure is transmitted into the micro-vascular networks, which can lead to pressure imbalance and vessel destruction (Feihl et al., 2009). Therefore, the BCW and BDW in the WI analysis can help to predict the loss of small vessels in the vascular beds (territorial ischemia). Although the mechanism of rarefaction-induced hypertension is complex (Feihl et al., 2009), our model analysis may provide a quantitative link between vasculature reduction and WI.
Finally, we study the effects of change in geometry in our model. Our results show that the FDW is very sensitive to geometry variation, which can have a 3-fold increase if the size is reduced by 10%, and 5-fold if the size is increased by 10%. This suggests that the FDW can be a good indicator of size change. However, we should also treat it with caution, as small errors in geometry measurements may induce large errors in the FDW. The other quantities, however, are less sensitive to size perturbations.
We now mention the limitations of our approach. The 1D model is based on a simplified tapered geometry for a 3D in vivo coronary tree. The complicated biophysical interaction between the ventricle and coronary flow is modeled using a dynamic space-dependent feedback pressure on the large vessels. In addition, we have not considered dynamic feedback pressure within the structured-tree model, which will change the impedance of the terminal vessels. Moreover, the venous pressure and flow were neglected in our model.
In summary, a 1D model based on the CT scans of coronary circulation is developed. The model employs a structured-tree model for the coronary vascular beds and includes combined feedback pressure resulting from contraction of the heart wall as an additional term in the boundary condition at the linking terminal arteries to the vascular beds. Our model agrees well with previous published studies and data. In addition, we use WI analysis to quantify flow and pressure waves in the coronary arteries. Pathological conditions such as stiffened coronary arteries and vascular rarefaction, as well as changes in geometry and the ventricular-pressure feedback ratio, are also studied.
Data Availability
The dataset for this manuscript is not publicly available because it is limited to clinical and research use in The Second Affiliated Hospital of Nanjing Medical University. Requests to access the datasets should be directed to ZD.
Author Contributions
ZD conceptualized and designed the study. The need for a structured-tree model for the coronary vasculature was identified by XL and NH. NH provided the code that was modified by WC and HG. XY helped with measuring the clinical experiment and collecting data. XL and NH helped with the analysis and interpretation of the data, critically revising the manuscript, and adding important intellectual content. All authors gave approval for the final version of this manuscript to be published and agreed to be accountable for all aspects of the work.
Funding
This work is supported by the Program for Changjiang Scholars and Innovative Research Team in University (Grant No. IRT_16R07), the National Natural Science Foundation of China (Grant Numbers: 11321062), by China Postdoctoral Science Foundation Funded Project (2017M623296XB), and Guangxi Young and Middle-Aged Teachers' Research Fundamental Ability Enhancement Project (2019KY0098), and the UK Engineering and Physical Sciences Research Council (Grant Number EP/N014642/1).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
The authors thank the patient and hospital colleagues for their assistance.
References
Alastruey, J., Parker, K., Peiro, J., and Sherwin, S. (2008). Lumped parameter outflow models for 1-D blood flow simulations: effect on pulse waves and parameter estimation. Commun. Comput. Phys. 4, 317–336.
Algranati, D., Kassab, G. S., and Lanir, Y. (2010). Mechanisms of myocardium-coronary vessel interaction. Am. J. Physiol. Heart Circ. Physiol. 298:H861. doi: 10.1152/ajpheart.00925.2009
Avanzolini, G., Barbini, P., Cappello, A., and Cevenini, G. (1988). CADCS simulation of the closed-loop cardiovascular system. Int. J. BioMed. Comput. 22, 39–49.
Chen, W. W., Gao, H., Luo, X. Y., and Hill, N. A. (2016). Study of cardiovascular function using a coupled left ventricle and systemic circulation model. J. Biomechan. 49, 2445–2454. doi: 10.1016/j.jbiomech.2016.03.009
Chilian, W. M., and Marcus, M. L. (1982). Phasic coronary blood flow velocity in intramural and epicardial coronary arteries. Circ. Res. 50, 775–781. doi: 10.1161/01.RES.50.6.775
Chiribiri, A., Hautvast, G. L. T. F., Lockie, T., Schuster, A., Bigalke, B., Olivotti, L., et al. (2013). Assessment of coronary artery stenosis severity and location: quantitative analysis of transmural perfusion gradients by high-resolution MRI versus FFR. JACC Cardiovasc. Imaging 6, 600–609. doi: 10.1016/j.jcmg.2012.09.019
Davies, J. E., Whinnett, Z. I., Francis, D. P., Manisty, C. H., Aguado-Sierra, J., Willson, K., et al. (2006). Evidence of a dominant backward-propagating “suction” wave responsible for diastolic coronary filling in humans, attenuated in left ventricular hypertrophy. Circulation 113, 1768–1778. doi: 10.1161/CIRCULATIONAHA.105.603050
Duanmu, Z., Yin, M., Fan, X., Yang, X., and Luo, X. (2018). A patient-specific lumped-parameter model of coronary circulation. Sci. Rep. 8:874. doi: 10.1038/s41598-018-19164-w
Duncker, D. J., and Merkus, D. (2004). Regulation of coronary blood flow. Effect Coron. Art. Stenosis Arch. Mal. Coeur. Vaiss. 97, 1244–1250.
Feihl, F., Liaudet, L., and Waeber, B. (2009). The macrocirculation and microcirculation of hypertension. Curr. Hypert.Rep. 11:182. doi: 10.1007/s11906-009-0033-6
Hadjiloizou, N., Davies, J. E., Malik, I. S., Aguado-Sierra, J., Willson, K., Foale, R. A., et al. (2008). Differences in cardiac microcirculatory wave patterns between the proximal left mainstem and proximal right coronary artery. Am. J. Physiol. Heart Circ. Physiol. 295, H1198–H1205. doi: 10.1152/ajpheart.00510.2008
Hall, J. E. (2015). Guyton and Hall Textbook of Medical Physiology e-Book. Oxford, MI: Elsevier Health Sciences.
Hoffman, J., and Spaan, J. (1990). Pressure-flow relations in coronary circulation. Physiol. Rev. 70, 331–390.
Huo, Y., and Kassab, G. S. (2007). A hybrid one-dimensional/womersley model of pulsatile blood flow in the entire coronary arterial tree. Am. J. Physiol. Heart Circ. Physiol. 292, H2623–33. doi: 10.1152/ajpheart.00987.2006
Karimi, A., Navidbakhsh, M., Faghihi, S., Shojaei, A., and Hassani, K. (2013). A finite element investigation on plaque vulnerability in realistic healthy and atherosclerotic human coronary arteries. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 227, 148–161. doi: 10.1177/0954411912461239
Kim, H. J., Vignon-Clementel, I., Figueroa, C., Jansen, K., and Taylor, C. (2010). Developing computational methods for three-dimensional finite element simulations of coronary blood flow. Finite Elem. Analy. Design 46, 514–525. doi: 10.1016/j.finel.2010.01.007
Kyriacou, A., Whinnett, Z. I., Sen, S., Pabari, P. A., Wright, I., Cornelussen, R., et al. (2012). Improvement in coronary blood flow velocity with acute biventricular pacing is predominantly due to an increase in a diastolic backward-travelling decompression (suction) wave. Circulation 126, 1334–44. doi: 10.1161/CIRCULATIONAHA.111.075606
Lax, P., and Wendroff, B. (1960). CADCS systems of conservation laws. Commun. Pure Appl. Math. 13, 217–237.
Lee, J., Nordsletten, D., Cookson, A., Rivolo, S., and Smith, N. (2016). In silico coronary wave intensity analysis: application of an integrated one-dimensional and poromechanical model of cardiac perfusion. Biomechan. Model. Mechan. 15, 1535–1555. doi: 10.1007/s10237-016-0782-5
Mantero, S., Pietrabissa, R., and Fumero, R. (1992). The coronary bed and its role in the cardiovascular system: a review and an introductory single-branch model. J. Biomed. Eng. 14, 109–116.
Mittal, N., Zhou, Y., Ung, S., Linares, C., Molloi, S., and Kassab, G. S. (2005). A computer reconstruction of the entire coronary arterial tree based on detailed morphometric data. Ann. Biomed. Eng. 33, 1015–1026. doi: 10.1007/s10439-005-5758-z
Mynard, J. P. (2011). Computer Modelling and Wave Intensity Analysis of Perinatal Cardiovascular Function and Dysfunction. Ph.D. thesis, Melbourne.
Mynard, J. P., Penny, D. J., and Smolich, J. J. (2014). Scalability and in vivo validation of a multiscale numerical model of the left coronary circulation. Am. J. Physiol. Heart. Circ. Physiol. 306, H517–H528. doi: 10.1152/ajpheart.00603.2013
Mynard, J. P., and Smolich, J. J. (2015). One-dimensional haemodynamic modeling and wave dynamics in the entire adult circulation. Ann. Biomed. Eng. 43, 1443–1460. doi: 10.1007/s10439-015-1313-8
Mynard, J. P., and Smolich, J. J. (2016). Influence of anatomical dominance and hypertension on coronary conduit arterial and microcirculatory flow patterns: a multiscale modeling study. Am. J. Physiol. Heart Circ. Physiol. 311:H11. doi: 10.1152/ajpheart.00997.2015
Olufsen, M. S. (1999). Structured tree outflow condition for blood flow in larger systemic arteries. Am. J. Physiol. Heart Circ. Physiol. 276, H257–H268.
Olufsen, M. S., Hill, N. A., Vaughan, G. D., Sainsbury, C., and Johnson, M. (2012). Rarefaction and blood pressure in systemic and pulmonary arteries. J. Fluid Mechan. 705, 280–305. doi: 10.1017/jfm.2012.220
Olufsen, M. S., Peskin, C. S., Kim, W. Y., Pedersen, E. M., Nadim, A., and Larsen, J. (2000). Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions. Ann. Biomed. Eng. 28, 1281–1299. doi: 10.1114/1.1326031
Parker, K. H. (2009). An introduction to wave intensity analysis. Med. Biol. Eng. Comput. 47:175. doi: 10.1007/s11517-009-0439-y
Pietrabissa, R., Mantero, S., Marotta, T., and Menicanti, L. (1996). A lumped parameter model to evaluate the fluid dynamics of different coronary bypasses. Med. Eng. Phys. 18, 477–484.
Pijls, N. H., and De Bruyne, B. (2013). Coronary Pressure, Vol 195. Dordrecht: Springer Science & Business Media.
Qureshi, M. U., and Hill, N. A. (2015). A computational study of pressure wave reflections in the pulmonary arteries. J. Math. Biol. 71, 1525–1549. doi: 10.1007/s00285-015-0867-2
Raff, G. L., Abidov, A., Achenbach, S., Berman, D. S., Boxt, L. M., Budoff, M. J., et al. (2009). SCCT guidelines for the interpretation and reporting of coronary computed tomographic angiography. J. Cardiovasc. Comput. Tomogra. 3, 122–136. doi: 10.1016/j.jcct.2009.01.001
Sen, S. (2013). Assessment of Intra-Coronary Pressure and Flow Velocity Relations Distal to Coronary Stenoses to Derive a Novel Index of Stenosis Severity. Ph.D. thesis, Imperial College London, London.
Sen, S., Petraco, R., Mayet, J., and Davies, J. (2014). Wave intensity analysis in the human coronary circulation in health and disease. Curr. Cardiol. Rev. 10, 17–23. doi: 10.2174/1573403X10999140226121300
Siebes, M., Kolyva, C., Verhoeff, B.-J., Piek, J. J., and Spaan, J. A. (2009). Potential and limitations of wave intensity analysis in coronary arteries. Med. Biol. Eng. Comput. 47, 233–239. doi: 10.1007/s11517-009-0448-x
Spaan, J. A. (2012). Coronary Blood Flow: Mechanics, Distribution, and Control, vol 124. Amsterdam: Springer Science & Business Media.
Spaan, J. A., Cornelissen, A. J., Chan, C., Dankelman, J., and Yin, F. C. (2000). Dynamics of flow, resistance, and intramural vascular volume in canine coronary circulation. Am. J. Physiol. Heart. Circ. Physiol. 278, H383–H403. doi: 10.1152/ajpheart.2000.278.2.H383
Stergiopulos, N., Young, D., and Rogge, T. (1992). Computer simulation of arterial flow with applications to arterial and aortic stenoses. J. Biomechan. 25, 1477–1488.
Sun, Y.-H., Anderson, T. J., Parker, K. H., and Tyberg, J. V. (2000). Wave-intensity analysis: a new approach to coronary hemodynamics. J. Appl. Physiol. 89, 1636–1644. doi: 10.1152/jappl.2000.89.4.1636
Toyota, E., Ogasawara, Y., Hiramatsu, O., Tachibana, H., Kajiya, F., Yamamori, S., et al. (2005). Dynamics of flow velocities in endocardial and epicardial coronary arterioles. Am. J. Physiol. Heart Circ. Physiol. 57:H1598. doi: 10.1152/ajpheart.01103.2003
Uus, A., Liatsis, P., Jawaid, M. M., Rajani, R., and Benderskaya, E. (2015). “Assessment of stenosis introduced flow resistance in CCTA-reconstructed coronary arteries,” in Signals and Image Processing (IWSSIP), 2015 International Conference on Systems (London: IEEE), 313–320. doi: 10.1109/IWSSIP.2015.7314238
Van der Horst, A., Boogaard, F. L., van't Veer, M., Rutten, M., Pijls, N. H., and van de Vosse, F. N. (2013). Towards patient-specific modeling of coronary hemodynamics in healthy and diseased state. Comput. Mathemat. Methods Med. 2013:393792. doi: 10.1155/2013/393792
Vignon-Clementel, I. E., Figueroa, C. A., Jansen, K. E., and Taylor, C. A. (2006). Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Comput. Methods Appl. Mechan. Eng. 195, 3776–3796. doi: 10.1016/j.cma.2005.04.014
Weber, T., Auer, J., O'Rourke, M. F., Kvas, E., Lassnig, E., Berent, R., et al. (2004). Arterial stiffness, wave reflections, and the risk of coronary artery disease. Circulation 109, 184–189. doi: 10.1161/01.CIR.0000105767.94169.E3
Keywords: hemodynamics, coronary arteries, blood flow, pressure, wave intensity, vascular rarefaction, structured tree model
Citation: Duanmu Z, Chen W, Gao H, Yang X, Luo X and Hill NA (2019) A One-Dimensional Hemodynamic Model of the Coronary Arterial Tree. Front. Physiol. 10:853. doi: 10.3389/fphys.2019.00853
Received: 11 February 2019; Accepted: 20 June 2019;
Published: 09 July 2019.
Edited by:
Pasquale Pagliaro, University of Turin, ItalyReviewed by:
Alun Hughes, University College London, United KingdomTimothy W. Secomb, University of Arizona, United States
Copyright © 2019 Duanmu, Chen, Gao, Yang, Luo and Hill. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zheng Duanmu, ariesdmz@sina.com; Weiwei Chen, chenww2010@hotmail.com