Coronary blood flow is regulated to match the oxygen demand of myocytes in the heart wall. Flow regulation is essential to meet the wide range of cardiac workload. The blood flows through a complex coronary vasculature of elastic vessels having nonlinear wall properties, under transmural heterogeneous myocardial extravascular loading. To date, there is no fully integrative flow analysis that incorporates global and local passive and flow control determinants. Here, we provide an integrative model of coronary flow regulation that considers the realistic asymmetric morphology of the coronary network, the dynamic myocardial loading on the vessels embedded in it, and the combined effects of local myogenic effect, local shear regulation, and conducted metabolic control driven by venous O2 saturation level. The model predicts autoregulation (approximately constant flow over a wide range of coronary perfusion pressures), reduced heterogeneity of regulated flow, and presence of flow reserve, in agreement with experimental observations. Furthermore, the model shows that the metabolic and myogenic regulations play a primary role, whereas shear has a secondary one. Regulation was found to have a significant effect on the flow except under extreme (high and low) inlet pressures and metabolic demand. Novel outcomes of the model are that cyclic myocardial loading on coronary vessels enhances the coronary flow reserve except under low inlet perfusion pressure, increases the pressure range of effective autoregulation, and reduces the network flow in the absence of metabolic regulation. Collectively, these findings demonstrate the utility of the present biophysical model, which can be used to unravel the underlying mechanisms of coronary physiopathology.
Coronary blood flow is regulated to match myocardial metabolic demand under a broad range of physical activity (Duncker and Bache, 2008). Flow regulation during high metabolic demand (e.g., exercise, stress, etc.) is achieved via a feed-forward mechanism mediated by sympathetic β-adrenoceptor vasodilation, in combination with local feedback mechanisms (Tune et al., 2002). Sympathetic activation during exercise was found to stimulate cardiac myocytes (resulting in higher heart rate [HR] and contractile force) and to induce coronary vasodilation (Gorman et al., 2000b, 2010). This feed-forward mechanism was estimated to account for ∼25% of coronary exercise hyperemia (Gorman et al., 2000b). Importantly, though, cardiac denervation or pharmacological inhibition of autonomic control does not result in ischemia during exercise, indicating that other vasodilator mechanisms act to regulate blood flow (Duncker and Bache, 2008).
Flow control is achieved by diameter regulation of resistance vessels that range from arterioles to small arteries (<400 µm in diameter) obtained by vasoactivity of vascular smooth muscle cells (VSMCs; Goodwill et al., 2017). Three major mechanisms of VSMC regulation have been proposed: (1) myogenic control (Johnson, 1980; McHale et al., 1987; Kuo et al., 1990a; Miller et al., 1997), whereby vessels vasoreact in direct response to trans-luminal pressure; (2) shear-induced control (Holtz et al., 1984; Kuo et al., 1990b, 1995; Jones et al., 1995a), whereby tone is regulated by nitric oxide (NO) synthesized by vascular endothelial cells (ECs) in response to wall shear stress; and (3) metabolic regulation that affects small arterioles (Feigl, 1983; Kanatsuka et al., 1989; Jones et al., 1995b) via metabolic endogenously produced vasoactive agents. Previous studies have shown that vascular responsiveness to these three regulation mechanisms is heterogeneous where (1) the adenosine (a putative vasodilator) effect is significant in downstream microvessels (diameters <150 μm) and decreases with increase in diameter (Kanatsuka et al., 1989; Jones et al., 1995b; Kuo et al., 1995), (2) shear-induced dilation is more prominent in upstream larger microvessels (Kuo et al., 1995), and (3) myogenic control is negligible in small microvessels because of their relatively low intravascular pressures (Chilian et al., 1989; Tiefenbacher and Chilian, 1998).
Coronary blood flow control raises important questions regarding the impact of each regulatory mechanism and the interactions thereof as well as the role of extravascular loading by the surrounding myocardium in flow regulation (Muller et al., 1996). Clinically, the circumstances that lead to failure of flow regulation remain unclear, e.g., failure of the system to maximally vasodilate even under profound ischemia (Hoffman and Spaan, 1990). Such questions cannot be experimentally studied in vivo because of the difficulties in separating the vascular response between shear and pressure, and in measuring in vivo flow and pressure in the deeper layers of the beating heart. Although model simulation is a powerful approach to study the coronary circulation, there are several major challenges. The molecular and transport aspects of the flow regulation mechanisms are still inadequately understood (Feigl, 1983; Chen and Popel, 2006; Dash and Bassingthwaighte, 2006). Previous studies adopted empirical single-vessel models based on isolated vessel studies (Liao and Kuo, 1997). The heterogeneities in coronary vasculature (VanBavel and Spaan, 1992; Kassab et al., 1993), and in intramyocardial pressure and vascular responsiveness to regulatory mechanisms, require a simulation that considers regulation mechanisms in each arterial segment as part of an integrated dynamic flow analysis in a realistic coronary vasculature.
Previous studies of coronary flow regulation (Liao and Kuo, 1997; Cornelissen et al., 2002; Carlson et al., 2008) did not incorporate important network characteristics. First, the network structure considered was assumed to be symmetric, and orders were assigned to vessel groups such as “small arterioles,” “large arterioles,” and “small arteries.” In addition, these models considered isolated networks having no interaction with the surrounding myocardium. Finally, flow analysis was steady-state without considering the dynamic effect of coronary flow. It has been conjectured (DeFily and Chilian, 1995) that in coronary circulation, the transmural variations of myocardium–vessel interaction (MVI) coupled with the network asymmetric structure, and the interactions between the various heterogeneous flow regulation mechanisms, are likely to result in a highly complex flow that may allow for adequate myocardial perfusion even under extreme metabolic conditions.
The objective of this study is to integrate local feedback flow regulation models using a structure-based multiscale framework for the analysis of dynamic flow regulation in a realistic coronary network subject to heterogeneous transmural variation of the extravascular loading by the contracting myocardium, and incorporating the longitudinal heterogeneity of passive and active vascular mechanical properties. The model can elucidate the impact of each regulation mechanism, the nature of the interaction between them, and the interaction of flow regulation with the extravascular myocardial contraction. The model predictions are compared with published data.
Materials and methods
The network structure
A morphological-based (Kassab et al., 1993) topological reconstruction of left circumflex coronary arterial tree (Kaimovitz et al., 2005) was applied for flow analysis. To reduce computational effort, a subendocardial subtree containing 400 vessels from order 6 to order 0 was pruned from the left circumflex coronary tree. The subtree contained 195 bifurcations, 3 trifurcations, and 79 terminal order 1 vessels. Diameters were randomly assigned to the tree vessel segments based on global statistics of the measured diameters under a pressure of 80 mmHg (Kassab et al., 1993). Preliminary flow analysis revealed that this method of diameter assignment resulted in excessive heterogeneity in both perfusion and wall shear stress. Hence, a recently developed iterative diameter reassignment was implemented to reduce flow heterogeneity to be on par with published data (Austin et al., 1994; Kassab and Fung, 1995). Reassignment is based on the assumption that the flow in each upstream network vessel should be proportional to the number of terminal vessels it feeds. Iteration was terminated when perfusion heterogeneity (indexed by the coefficient of variation [CV] = SD/mean) was reduced to within the range of published data.
To study the effects of transmural location on the network flow, the same subendocardial subtree was “implanted” in the subepicardium layer. The network flow analysis was then subjected to the subepicardium boundary conditions (i.e., inlet and outlet pressures and extravascular loading by the MVI).
Network flow analysis
Details of the flow analysis are provided in Appendix A. Briefly, the flow in each elastic vessel was modeled by a validated (Jacobs et al., 2008) three-element Windkessel (Fig. 1, A and B) made of two nonlinear resistors in series and one parallel capacitor. Flow resistance was governed by Poiseuille’s equation and vessel capacitance by the vessel pressure–diameter relationship (PDR). Network flow was solved by imposing flow continuity at each vessel midpoint and at each network junction i.e., the net flow at each junction vanishes to conserve mass. This formulation resulted in a system of ordinary differential equations (ODEs) that were numerically solved (Appendix A). The network dynamic flow was analyzed subject to boundary conditions of network dynamic inlet and outlet pressures, and of the extravascular loading by the cyclic myocardium contractions.
Under the myogenic and shear regulation, network flow was iteratively solved by adjusting each vessel diameter according to the local pressure and flow within it. For cases of metabolic regulation aimed to achieve a set level of perfusion (matching the metabolic demand), flow was solved by optimizing the distribution of metabolic activation which produced terminal flows close to this perfusion level.
Vascular mechanical properties
The passive and active vessel PDR depends both on the vessel order and diameter (Liao and Kuo, 1997). Published data are available for selected diameters. In the current work, the passive and active properties of all pertinent vessel sizes were deduced from the available data and extrapolated where needed.
The passive vessel properties
The passive vessel response consists of two contributions: the intrinsic vessel passive mechanics (represented by its PDR) and the tethering effect of the myocardial tissue in which the vessel is embedded. The PDR of isolated in vitro passive vessels under positive transvascular pressure showed the radius to be a sigmoidal function of the dynamic transvascular pressure, ΔP (Young et al., 2012) as
where Ap and Bp are the asymptotical highest and lowest radii, respectively (under the highest positive and lowest negative transvascular pressure, respectively), is the transvascular pressure corresponding to the mean of radii Ap and Bp, and Cp is the passive response bandwidth. The tethering model is described in Appendix C (Eqs. C1, C2, C3, and C4) and the data on the PDR constants for orders 1–6 vessels are listed in Table A2.
Active vessel properties
The active wall properties are determined by VSMCs, which are affected predominantly by myogenic (pressure), flow (shear), and metabolic control mechanisms. Published information on vascular response to these mechanisms relates primarily to the steady-state response (Liao and Kuo, 1997). There is paucity of information on the dynamics of coronary control, and much of the transient response relates to changes in HR including the combined effects of multiple flow regulation mechanisms (Dankelman et al., 1992). The analysis of flow regulation, however, requires detailed knowledge of the dynamics of each flow control mechanism, which is currently unavailable. Here, we focused on flow features that have been experimentally studied under steady-state responses including the functional role of each control mechanism and of their interactions. Although the present study focused on the steady-state active response under constant HR, the flow and mechanics of the vascular system and its interaction with the myocardium were dynamically analyzed under these quasi-static control conditions (see The flow regulation time constant section below).
The myogenic regulation: The myogenic regulation results from VSMC contraction in response to local wall stress as determined by the transvascular pressure. The myogenic diameter reduction is expressed by a sigmoidal function of the time-averaged transvascular pressure (Young et al., 2012) as
where ρm is the myogenic response amplitude, is the transvascular pressure under which the myogenic radius change is highest, Cm is the myogenic response bandwidth, and m is a shape factor. The procedure of assigning parameters for Eq. 2 to the different vessel orders and their values are given in Appendix F and Table A3, respectively.
The regulated vessel radius under quasi-static loading was taken to be a function of the total activation level, A, and of the mean transvascular pressure, (Liao and Kuo, 1997). The maximum myogenic reduction in radius ΔRm under full activation (A = 1), is attenuated with the total activation (A < 1; Liao and Kuo, 1997) to yield the regulated radius as
where the total activation (A) caused by myogenic, shear stress and metabolites is given by
Expressions for Fτ and Fmeta are listed below (Eqs. 5 and 6). The product form in Eq. 4 between the metabolic and flow regulations pertains to their respective residual activities (1 − Fτ) and (1 − Fmeta). Importantly, this form is mathematically identical to the physiologically based additive model proposed and experimentally validated by Liao and Kuo (1997).
Flow (shear) regulation: The shear regulation induces relaxation of the myogenic contracted vascular wall, which is mediated by NO production by the ECs in response to local wall shear stress. The shear fractional deactivation is taken to be dependent on the mean shear stress, and is expressed by Liao and Kuo (1997):
where Kτ is the wall shear stress constant and Fτmax is the maximum deactivation caused by wall shear stress. The procedure of assigning parameters for Eq. 5 to the different vessels’ orders is summarized in Appendix F, and the parameters are listed in Table A4.
Liao and Kuo (1997) pointed out that in modeling in vivo network flow, the values of their in vitro measured Kτ are too low because of the presence of hemoglobin, which binds to NO and hence decreases the in vivo sensitivity to shear. In their flow analysis in an idealized symmetric network without MVI, they increased Kτ by a factor of 150, which allows the vessels to respond to shear stress under physiological conditions. In the present network simulations, it was found that a factor of 15 is sufficient (Table A4).
The metabolic regulation: Based on the evidence presented in Appendix E, the mechanism of metabolic regulation considered here is the conducted response (CR). In view of the coronary vascular network structure, the model of CR integrates the signaling effects along the various network pathways, from the precapillary arterioles to each specific upstream vessel. Hence, the metabolic activation in an upstream vessel i is expressed by
where ui is the number of terminal vessels fed by the upstream vessel i, is the strength of the response in an ith vessel conducted from its jth terminal vessel, and L0 is the decay characteristic length. The metabolic signal in the jth terminal order 1 vessel is a direct function of the local oxygen supply/demand imbalance. The CR was found to decay exponentially with the path length, , toward the upstream vessels (Delashaw and Duling, 1991; Xia and Duling, 1995; Arciero et al., 2008; Goldman et al., 2012) with a characteristic length, L0 which determines the rate of decay of the metabolic activation with path length (Eq. 6). The total metabolic activation, , in an ith vessel was taken to be the mean from all jth terminal vessels fed by that vessel. A reference value of L0 = 1 mm was selected for the decay characteristic length. This value lies within the measured range of different vascular beds (0.15–2.5 mm; Hald et al., 2012).
The metabolic signal in each j terminal order 1 vessels, , depends on the local demand/supply imbalance. In control theory, the control signal is the system desired output, which in the coronary circulation is the requisite terminal arteriole perfusion (qtarget) that balances the metabolic O2 demand. This is irrespective of the metabolites or mechanisms involved. Hence, although the terminal arteriole flow is not physiologically a sensed signal, its requisite level represents the metabolic demand regardless of the involved metabolite pathways. Our choice of a metabolic signal is supported by the findings that coronary flow correlates well with an increase in the coronary venous ATP concentration, and the latter correlates with the decline of venous PO2 (Farias et al., 2005). In the network flow analysis, the terminal vessel metabolic signal was either set to be constant in all terminal vessels or optimized in each terminal vessel to provide a set level of terminal flow.
The flow regulation time constant: The time constant of coronary vessel response to changes in pressure and flow (∼1.5-fold of t50, the time required to establish half of the complete response) was found to be in the range of 15 s to 2 min (Mosher et al., 1964; Hoffman and Spaan, 1990; Dankelman et al., 1992; Tsoukias et al., 2004). This response time constant is significantly higher than the cardiac period (∼1 s). The stabilized system response can thus be considered as the time-average over a cardiac cycle. Hence, the levels of regulated vessel radius (Eq. G1), of the active tension (Eq. G4), and of the active stiffness (Eq. D5) were formulated as functions of the time-averaged transvascular pressure.
The boundary conditions: Algranati et al., 2010 showed that in addition to the inlet and outlet pressures, the coronary flow is determined by the MVI, which consists of the combined effect of the intramyocardial fluid pressure (IMP) and the shortening-induced intramyocyte pressure PSIP. IMP varies with the myocardial relative depth (MRD) from the left ventricle (LV) pressure at the endocardium to zero at the epicardium. Waveforms of the inlet pressure, Pin(t), outlet pressure, Pout(t), LV pressure, PLV(t), and intramyocyte pressure, PSIP(t), are input signals to the flow analysis (Fig. A1). The Pout(t) signal was interpolated for different transmural locations based on predictions from simulation of the unregulated flow in an entire coronary network that included arterial and venous trees and four identical representative capillary networks, at MRDs of 0.125, 0.375, 0.625, and 0.875 (Algranati et al., 2010). PLV(t) waveform was taken from predictions based on a preliminary study of a distributive LV mechanical model under resting HR (Kiyooka et al., 2005) of 75 BPM. Several considerations guided the choice of the Pin(t) signal for the subendocardial 400-vessel network. The first is the pressure drop from the aorta to the trunk vessel (order 6) of the subtree. On the other hand, there is a pressure increase caused by the added intramyocyte pressure, PSIP(t), which develops during contraction (Rabbany et al., 1989). Finally, Pin(t) must provide for sufficient flow perfusion in the terminal order 1 vessels in the range of measured flow of 0.4–2.0 × 10−3 mm3/s in systole and diastole (Tillmanns et al., 1974). Based on these considerations, Pin(t) was chosen to be 122/90 mmHg (with mean = 100) in systole/diastole, and the signal shape was adopted from (Algranati et al., 2010). Pout was assigned for each terminal vessel to be between the previously predicted subepicardium and subendocardium signals and (Fig. A1; Algranati et al., 2010), depending on the transmural location of the vessel. In the absence of data on Pout changes under higher metabolic demands (i.e., higher qtarget), the level for each vessel was kept the same under changes of qtarget. The tissue pressure PT(t) was derived based on the previous analysis of unregulated coronary flow (Algranati et al., 2010).
Model comparison with flow characteristics
Because detailed quantitative validation data are presently not available, predictions were compared against global response patterns. To that end, four sets of simulations were performed under the following conditions: (1) full myogenic activation; (2) myogenic and shear activation with no metabolic activation; (3) myogenic, shear, and full metabolic activation of terminal order 1 vessels; and (4) myogenic, shear, and optimized metabolic activation in terminal order 1 vessels aimed to achieve target flow levels in the terminal arterioles. Case 1 is an arrested heart with the sole effect of pressure on flow control. The shear and metabolic signals were deactivated by setting all Fτ and all Fmeta in Eq. 4 to zero. Case 2 is also an arrested heart with the effects of pressure and shear on flow control. The metabolic signal was deactivated by setting all Fmeta in Eq. 4 to zero. Case 3 represents the highest flow capability of the system in a beating heart. In case 4, the optimal metabolic signals, , of all terminal order 1 vessels were targeted to obtain the least deviation of terminal perfusion level, , from the set target flow level (qtarget). The details of the flow simulation are presented in Appendix A. A listing of the reference case run and model parameters is given in Table 1.
The following model predictions were compared with observations: dispersion of the perfusion, transmural perfusion heterogeneity, the magnitude of the metabolic flow reserve (MFR), and autoregulatory response of the system flow.
Perfusion dispersion: Coronary perfusion in the passive (unregulated) state has been shown to be highly heterogeneous (Austin et al., 1990; Huo et al., 2009). Regulation was found to decrease flow dispersion (Austin et al., 1994). The predicted flow level and its dispersion (CV) were evaluated here under various regulation mechanisms with full (Fmterm = 1) and no (Fmterm = 0) metabolic activation (Fig. 2 A and Table 2), and under optimized metabolic activation for several levels of metabolic demand (i.e., target terminal vessel flow, qtarget) and a few levels of input perfusion pressures (Fig. 2 B and Table 3).
Transmural perfusion heterogeneity: The effect of network transmural location on perfusion and its dispersion was analyzed as (1) terminal arteriole flow qterm under various regulation mechanisms and reference perfusion pressure and (2) optimized flow under a few levels of metabolic demand and perfusion pressures (Tables 2 and 3). Transmural heterogeneity was also studied in terms of the MFR, the network total flow under different regulation mechanisms, and autoregulation of the network flow.
Metabolic flow reserve: MFR is the flow reserve that results from flow regulation alone, without the effects of HR and activation waveform. It was estimated (Table 4) from the total flow in the active vessel network under myogenic, shear, and metabolic activations. MFR was estimated (Fig. 3) as a function of the mean input pressure () in the range of 45–180 mmHg under the same waveform and outlet pressure signal ().
Flow autoregulation: Autoregulation in coronary circulation under constant metabolic demand produces approximately constant network flow under an increase in perfusion pressure. It was studied with the metabolic regulation optimized to provide four qtarget levels of 1.20, 1.50, 1.80, and 2.25 × 10−3 mm3/s and under increasing perfusion pressure. For comparison, the flow was analyzed also for the passive state, and under full (Fmterm = 1) and no (Fmterm = 0) metabolic activation (Fig. 3).
Effect of MVI: The vessel loading by the myocardium (MVI) has three passive components (Algranati et al., 2010; Young et al., 2012): (1) chamber derived intramyocardial tissue pressure, (2) intracellular pressure in the contracting myocytes, and (3) effect of vessel tethering to the surrounding myocardium. The effects of MVI on the flow were analyzed in terms of the MFR (Table 4), the significance of each regulation on the total network flow (Fig. 3), and flow autoregulation (Fig. 4), and its effect on the distribution of metabolic activation for different levels of metabolic demand and input perfusion pressure (Fig. 5).
Flow sensitivity to each regulation mechanism
Flow sensitivity to each regulation mechanism was studied in two different ways. First, it was deduced from the effect of each regulation mechanism on the relationship between the network flow and the input perfusion pressure , i.e., from the vertical shift of the curve under each mechanism. This was done both with and without the effect of myocardium extravascular loading (MVI), and for both the subendocardial and subepicardial networks (Fig. 3). Second, sensitivities to the myogenic and to the shear regulations were investigated by perturbation of their respective amplitudes (ρm and Fτmax) to both higher and lower levels in reference to the baseline (Table 5).
Online supplemental material
Fig. S1 depicts the tree structure. Table S1, included in a separate PDF, shows the vessel dimensions and the network connectivity.
The integrated myogenic and flow regulation mechanisms, combined with the conducted metabolic response, reproduce the experimentally observed flow features including significant flow reserve, autoregulation of flow under a range of perfusion pressures, and lower dispersion of regulated flow except under extreme levels of metabolic demand. The predictions are also in accordance with observations that the myogenic and metabolic flow regulations are dominant, with lesser effect of the shear regulation. These findings are consistent with observed features of the coronary system (e.g., Kanatsuka et al., 1989; Austin et al., 1994; Matsumoto and Kajiya, 2001).
CV in precapillary flow rates is a measure of flow dispersion (Tables 2 and 3). Under the reference inlet and outlet pressures (Table 1) and passive vessel conditions, and under full metabolic activation, CV was ∼25% and ∼29%, respectively (Table 2). Under just myogenic regulation it increased to ∼35% and 46% in the subendocardial and subepicardium networks, respectively. When the metabolic regulation was optimized to yield the target terminal flow, the flow dispersion reduced to between ∼10% and ∼30%, depending on perfusion pressure and target flow (Table 3). In general, the flow dispersion is high under extreme conditions where the metabolic demand (qtarget) is low and the perfusion pressure is high, and when the metabolic demand is high and perfusion pressure is low (Table 3). In these cases, the metabolic control is either shut down (all Fmterm = 0) or exhausted to its full capacity (all Fmterm = 1) and the flow CV is very high (see also Fig. 5).
Transmural perfusion heterogeneity
The analysis revealed transmural perfusion differences, as exhibited by the terminal flow rate and dispersion, by the MFR and the effects of MVI, and by the effects of metabolic demand (indexed by the target terminal flow) and perfusion pressure. The predicted flow rate is higher in the subepicardium under all regulation conditions (Table 2). These transmural flow differences reduce under optimized metabolic regulation aimed to achieve a set target flow (Table 3) to levels that are flow rate dependent. Under optimized metabolic activation (Table 3), the CV is ∼10% in both the subendocardial and subepicardium networks for the reference target terminal flow of 1.50 × 10−3 mm3/s. This transmural trend changes significantly, however, with the target flow and perfusion pressure. Under all target flow rates studied, the dispersions are lower in the subepicardium than in the subendocardium under low perfusion pressures, but are higher than in the subendocardium under high perfusion pressures (Table 3).
Metabolic flow reserve
MFR is variably affected by the input perfusion pressure depending on the network transmural location and MVI (Table 4). In the in vivo case (with MVI) under low perfusion pressure, it is higher in the subepicardium, whereas under medium and high , it is higher in the inner subendocardial layer. In the subepicardium, MFR had a similar trend as in the subendocardium under no MVI (Table 4). In the subepicardium, MFR decreases monotonously with increasing perfusion pressure from 3.25 at = 60 mmHg to 1.45 at = 180 mmHg. In the subendocardium, MFR is highest (3.19) under = 120 mmHg, and drops under both lower and higher perfusion pressures (Table 4).
Effect of regulation on the coronary flow
The total flow increases with the mean inlet pressure () under all flow regulation mechanisms, both with and without MVI (Table 3 and Fig. 3). Comparison of flow distributions in the precapillary (order 1) terminal vessels under four flow regulation conditions (Fig. 2 A) shows that under solely myogenic activation (when the vessel diameters are at their lowest values), the distribution shifts to the extreme left of low flow levels. Under myogenic and shear regulations (no metabolic regulation), the distribution shifts to slightly higher flow levels. Full metabolic activation shifts the flow distribution to the right (toward high flow levels), closer to their passive levels.
Under optimized metabolic activation with the reference level of perfusion pressure, the mean terminal flow becomes closer to the set target flow qtarget (Fig. 2 B) but deviates significantly from it under both lower and higher perfusion pressures (Table 3). As expected, comparison with the results in Table 2 indicates that a significantly lower flow CV is obtained under optimized metabolic activation (set level of qtarget) as compared with cases of either full or no metabolic activation.
Simulations show that autoregulation occurs under increasing , as evident from the flow plateau region in the flow-pressure curves (Fig. 4 A). The pressure range where autoregulation is effective is higher under higher target flows. In the absence of MVI (Fig. 4 B), the flows under all conditions increase as compared with the flow with MVI. Also, MVI increases the effective autoregulation pressure range for all target flows. In the subepicardium network with MVI, the flow-pressure curves (Fig. 4 C) under all conditions (including autoregulation) are similar to those of the subendocardial network without MVI (Fig. 4 B).
The distribution of the metabolic activation (Fmterm) under flow optimization and reference target terminal flow qtarget of 1.5 × 10−3 mm3/s depends on the input perfusion pressure. The distribution shifts to lower values under high (Fig. 5, E and F). On the other hand, the metabolic activation in the terminal vessels is at the maximum level (all Fmterm ≃ 1) when inlet pressure is low ( = 75 mmHg). The distribution moves toward the lowest possible level (all Fmterm ≃ 0) as increases and becomes zero in all vessels when is >135 mmHg (Fig. 5, F and K).
Effect of MVI
MVI enhances the effect of myogenic and metabolic regulation in the sense that it induces a higher reduction in the myogenic flow and increases metabolic flow recovery. But it has a minor effect on the shear regulation (Fig. 3). MVI has a significant and -dependent effect on the MFR in the subendocardial network (Table 4). Under low perfusion pressure (60–90 mmHg), MVI reduces the level of MFR. Under medium to high pressure (100–180 mmHg; Table 4), MVI significantly increases the MFR level. MVI has a much lower effect on the subepicardium network flow. This is evident from the closeness of the MFR levels in the subepicardium with MVI to that of the subendocardium without MVI except under very low perfusion pressures (Table 4).
Flow sensitivity to each regulation mechanism
The different control mechanisms have a variable effect on the local flow dispersion (Table 2). Flow CV is highest under just myogenic regulation, followed by myogenic + shear, myogenic + shear + full metabolic activation, and under no regulation (passive state).
The network flow was found to be highly sensitive to the myogenic regulation (Table 2 and Fig. 3), as evident from the significant reduction in flow levels as compared with the passive state, and from the large effect of the myogenic amplitude ρm (Eq. 2) on the flow (Table 5). The mean flow in terminal order 1 vessels decreased from 0.30 × 10−3 mm3/s to 0.08 × 10−3mm3/s when ρm was increased by 40% of the reference level and increased to 1.05 × 10−3mm3/s when ρm was decreased by 50% from its reference level. Flow dispersion is affected as well. The CV increased from 0.27 to 0.36 when ρm increased by 90%.
The predicted sensitivity to shear is lower than that for the myogenic regulation (Fig. 3). The mean terminal flow decreased from 0.63 × 10−3 to 0.44 × 10−3 mm3/s when Fτmax was reduced to 50% of the reference and increased to 0.83 when Fτmax increased by 40% of its reference level. The corresponding flow CV decreased slightly with increasing Fτmax, and vice versa (Table 5).
There is a significant sensitivity of the flow to the metabolic regulation as seen by the large shift of the network flow versus perfusion pressure curves to close to the passive curves under full metabolic activation (Fig. 3). The effect of optimized metabolic activation increases with increasing target flow, to an extent that depends on the perfusion pressure and on the network location and MVI (Fig. 4).
Order dependence of the metabolic diameter regulation
The model predicts that as the metabolic demand increases from a low level (qtarget = 1.2 × 10−3 mm3/s) to a high level (qtarget = 2.25 × 10−3 mm3/s) under perfusion pressure of 120 mmHg, the mean diameter increase is 15% for terminal order 1, 13% for order 2, 11% for order 3, 10% for order 4, and 9.5% for orders 5 and 6 vessels.
Novel model outcomes
Two new outcomes are gleaned from the model. First is failure of the regulation system to maximally vasodilate the vessels even under profound ischemia. Second, the extravascular loading by the contracting myocardium (MVI) enhances the flow reserve under high perfusion pressure and extends the range of perfusion pressure under which the flow is effectively autoregulated. These new findings could not have been captured previously. In particular, the second insight on the effect of MVI is counterintuitive. The rationale for both findings is presented in the Discussion.
The present biophysical platform provides a rational basis for the system characteristics and a coherent insight into the underlying mechanisms. It advances flow regulation in several respects as outlined below.
The validity of the present integrated simulation is based on the realism of its components: (1) the network structure was reconstructed (Kaimovitz et al., 2005) from measured morphometric data on coronary vasculature (Kassab et al., 1993; Kassab and Fung, 1994), (2) the effects of both the transmural extravascular myocardial loading and the vessel tethering to the myocardium are incorporated, as opposed to isolated, externally unloaded networks, (3) the models for the three regulation mechanisms were adopted from a detailed experimental study on coronary microvessels (Liao and Kuo, 1997), (4) the conducted nature of the metabolic regulation is a key element in the present study and was not considered in previous coronary models (Liao and Kuo, 1997; Cornelissen et al., 2002), and (5) the present study incorporates the dynamics of coronary flow as opposed to previous steady-state models (Liao and Kuo, 1997; Cornelissen et al., 2002). The dynamic aspects may not be obvious because results are presented in terms of time-averaged quantities, but they are indeed captured by the dynamic boundary conditions (Fig. A1) and the dynamic MVI. In addition, (6) the network dynamic flow analysis including the mechanisms for the MVI were previously validated against measured data (Algranati et al., 2010).
Model validation is made by comparison with several observed coronary flow characteristics. A comparison of the flow sensitivities to the three regulation mechanisms (Fig. 3 and Table 5) indicates that the myogenic and metabolic regulations have major effects on the flow, whereas the flow (shear) mechanism has a minor effect. These predictions are consistent with experimental observations (Feigl, 1983). Goodwill et al., 2017 presented a number of studies showing that NO-dependent responses occur primarily in upstream large arterioles and arteries (100–300 µm). The NO responses are largely absent in resistance vasculature (<100 µm diameter; Chilian et al., 1993), a diameter range that includes most of the resistance vessels in the network of the present study.
The important role of the metabolic regulation is clearly demonstrated in Figs. 3 and 4, where the metabolic mechanism is responsible for a number of essential coronary flow features. Of special note is the network’s ability to meet the flow levels required to balance the metabolic demand. In so doing, it also significantly reduces flow dispersion and transmural flow heterogeneity. Finally, the metabolic mechanism facilitates autoregulation under increasing perfusion pressure. The analysis has shown that the metabolic and myogenic regulations alone are sufficient to yield flow autoregulation of a similar pattern to that in Fig. 4 A (with a larger range of autoregulation pressure range) without the need for shear regulation.
The model prediction that under increased metabolic demand, the diameter of smaller arterial microvessels increases more than those of the higher order ones is in accordance with the experimental conclusion (Kanatsuka et al., 1989) from studies with dogs under increased oxygen consumption produced by rapid pacing.
Coronary flow reserve (CFR) is an important functional characteristic of the network, which has two major components. One derives from the metabolic regulation (the MFR). The other is a result of increased HR and altered waveform (reduced diastolic time fraction [DTF]) under high metabolic demand. The model predicts MFR of a factor of 2–3.5 under physiological levels of perfusion pressure (Fig. 3 and Table 4), whereas CFR can be as high as 5–6 in young healthy adults (Pitkänen et al., 1998). The combined effects of increased HR (factor of 2–3.5) and increased perfusion pressure are expected to meet the measured CFR.
Clinical observations indicate that even under profound ischemia, flow regulation fails to maximally dilate the vessels (Hoffman and Spaan, 1990; Sambuceti et al., 2001). This observation can be accounted for by recalling that even under full metabolic activation of the terminal vessels, the predicted flow is lower than the passive one (Fig. 3 and Table 2). Even if terminal vessels are maximally vasodilated by the metabolic signal, not all upstream vessels are maximally dilated because of the decaying nature with path length of the CR toward those upstream vessels.
Under normal physiological inlet perfusion pressure and metabolic demand, metabolic regulation yields the required metabolic needs (represented by qtarget) and enhanced perfusion homogeneity (Table 3, the reference case). Under extreme levels of demand or perfusion pressure (either very low or very high), however, there is a growing difference between predicted and target flows, as well as increased flow dispersion. The mechanistic basis is that under such extreme conditions, the metabolic regulation is either in full activation (all Fmterm = 1) when is extremely low (Fig. 5 A), or completely suppressed (all Fmterm = 0) when is extremely high (Fig. 5 F and Table 3). A similar metabolic limitation is observed in the autoregulation response (Fig. 4 A). Autoregulation has the highest effect in the physiological range of inlet pressure. As the inlet pressure falls outside of the physiological range, autoregulation cannot be maintained and the flow deviates from the metabolic demand (qtarget).
MVI affects the network pressure-flow relationship (Fig. 3) primarily under myogenic and myogenic + shear regulations but seems to have little effect on the passive and full metabolic flows. MVI also enhances the MFR (Table 4) for all cases except under low inlet perfusion pressures. Another important effect is on the flow autoregulation (Fig. 4), i.e., the presence of MVI significantly increases the pressure range of effective autoregulation. Under in vivo conditions (with MVI), the autoregulation pressure range is expected to be even wider than in Fig. 4 because the coronary perfusion pressure increases with increasing LV pressure. The latter is a key determinant of the magnitude of MVI extravascular loading.
The large spatial flow dispersions of ∼24% observed in the passive networks in both subendocardium and subepicardium compare well with the measured flow dispersion of ∼30% in a beating heart without tone (Austin et al., 1994). The predicted flow dispersion under physiological levels of perfusion pressures (∼10–15%) under autoregulation (with vascular tone) is slightly higher than 5–10% reported in one study (Matsumoto and Kajiya, 2001) but within the range of 14–19% reported in another (Austin et al., 1994). The differences between these experimental studies are likely the result of differences in methodology (double-tracer digital autoradiography vs. microsphere), species (rabbit vs. canine), and sample size (1 vs. 150 mg). In the present study, the estimated tissue volume perfused by one terminal arteriole is on the order of 1 mg. Extrapolation of flow dispersion from a large sample size (1 g; Matsumoto and Kajiya, 2001) to a small sample (1 mg) based on fractal scaling (Bassingthwaighte et al., 1989) yields a flow dispersion in the small sample that is significantly higher than that measured by Matsumoto and Kajiya, 2001. In addition, microsphere measurements are known to overestimate flow dispersion, and microspheres of 15 µm can interfere with the vessel wall and tone in the small arterioles. Higher flow dispersion can lead to local oxygen deficit and local hypoxic conditions in the myocardium. To summarize this point, the predicted flow dispersion of ∼10–15% under autoregulation seems reasonable, and in close agreement with data of Matsumoto (Matsumoto and Kajiya, 2001).
Transmural flow heterogeneity is found with higher precapillary network flow rates in the subepicardium as compared with the subendocardium (Table 2). This is attributed to the reduced subepicardium MVI allowing vessels to dilate to a greater degree than subendocardial vessels. During autoregulation, however, the network flow rate in the presence of MVI in the normal physiological state was found to be similar across transmural regions (Table 3 and Fig. 4, A and C). An additional important transmural difference is in the perfusion density (perfusion per tissue volume). This measure is affected by both the terminal flow and the vessel density in the myocardial tissue. Capillary density was found to be ∼30% higher in the subendocardium as compared with the subepicardium (Gerdes and Kasten, 1980; Breisch et al., 1986; Lee et al., 2009). Under in vivo autoregulation conditions, measured perfusion density was found to be ∼20% higher in the subendocardium than in the subepicardium (Feigl, 1983). The model predictions are in close agreement with these data showing that for three levels of metabolic demand (qtarget = 1.25, 1.5, 1.75 × 10−3 mm3/s), the endocardium/epicardium perfusion density ratios are 1.1, 1.26, and 1.26, respectively. In the passive (unregulated) state, the model predicts that the maximally achievable perfusion density under a perfusion pressure of 100 mmHg is slightly higher in the subendocardium than the subepicardium (1.37 vs. 1.24 ml/min/g), which follows a similar trend to the data of Fokkema et al. (2005) on the transmural distribution of the conductance.
The present study has several limitations. First, the networks analyzed are subendocardium and subepicardium small subtrees (400 vessels) compared with the full-depth coronary network. Analysis of the full network is preferred because such a network covers the entire transmural depth and can incorporate more accurately the transmural heterogeneity in MVI and vessel properties. The model boundary conditions may present another limitation. The pressure boundary conditions (Pin, Pout, and MVI) were adopted from mechanical analysis of a thick-wall cylindrical LV (Fig. A1), it is believed that these boundary conditions can realistically represent the equatorial region of the LV. Analysis of the entire coronary network requires boundary conditions that are either measured or estimated from a realistic model of a heart. However, in view of the local nature of the flow regulation and the realistic structure of the analyzed network (Kaimovitz et al., 2005), the approach of applying boundary conditions predicted by a thick-wall cylindrical model is reasonable. Another limitation is that the study pertains to a single HR and DTF representative of resting conditions. Flow regulation is essential under higher levels of physical activity characterized by higher HR and reduced DTF. In addtion, omission of the effect of collateral vessels (Schaper and Schaper, 1993) is another limitation of the present analysis. The flow regulation model has several limitations. First, the assumed metabolic feedback mechanism at the level of individual vessel has infinite gain within the dynamic diameter range of that vessel order. At the network level, however, the terminal flow is dispersed even when only a minority of the terminal order 1 and of the upstream vessels attain the limits of their dynamic diameter range. This is a result of the asymmetric branching pattern and transmural heterogeneity of the extravascular loading on the vessels. For example, under a perfusion pressure of 120 mmHg, the percent of terminal vessels with diameters that are different from their passive (highest) ones by >5% is 92%, 82%, 66%, and 32% under target terminal flows (qtarget) of 1.2, 1.5, 1.8, and 2.25 × 10−3 mm3/s, respectively. For the nonterminal vessels, the corresponding figures are 100%, 100%, 96%, and 82%. An additional limitation with the metabolic model is that it incorporates only upstream decaying CR. Because of a lack of sufficient experimental information, the model does not include the possibility of nondecaying CR described by Figueroa and Duling (2008).
A theoretical limitation of the model relates to the myogenic model (Eq. 2) that was adopted from the experimentally validated model of (Liao and Kuo, 1997) to ensure realism. That study on in vitro isolated vessels considered the intravascular pressure as the myogenic signal. Physiologically, it is the circumferential tension that is the major loading on the wall smooth muscle. It is hence likely the signal to which they respond by contracting, thus changing the vessel diameter. The wall tension depends on both the transvascular pressure (in a proportional manner) and on the vessel diameter. This theoretical limitation is circumvented here because in the data and model of (Liao and Kuo, 1997), dependence on diameter was implicitly incorporated via the diameter dependence of the myogenic parameters. Finally, the model validation is qualitative, which stems from a lack of sufficient data on model parameters in a specific heart with known structure, a requisite for quantitative predictions and subsequent comparison with data. Faced with this paucity of data, the approach was to compare observed response patterns to predictions of the model in which all the elements (such as stochastic asymmetric structure and regulation models of individual vessels) were taken from previous experimental studies, and augmented by interpolation when needed (Tables A2, A3, and A4).
Rationale of the new model outcomes
The failure of the regulation system to maximally vasodilate the vessels even under profound ischemia is a direct consequence of the conducted vasodilatory signaling, which implies that even under full metabolic activation of the terminal vessels, not all the network vessels are maximally dilated to their passive diameters.
The two novel and apparently counterintuitive model outcomes that the extravascular loading by the contracting myocardium (MVI) enhances the flow reserve under high perfusion pressure and extends the range of perfusion pressure under which the flow is effectively autoregulated are intriguing. The first, seen in Table 4, is unexpected in view of the compressive effect of MVI on the vessels. A closer inspection based on the conceptual framework of the present model reveals the rationale for this result. In the absence of MVI and under low metabolic demand, vessels are relatively dilated because of their high transvascular pressure, which provides little reserve for further dilation under increased metabolic demand. The resulting flow reserve (ratio of highest-to-sedentary perfusion levels) is thus low. With MVI, on the other hand, the transvascular pressure is reduced by the extravascular loading, and so are the vessel diameters. Hence, the vessels retain their capacity to significantly dilate under higher metabolic demand and hence increase the flow. The resulting flow reserve is thus increased.
The second prediction on the effect of MVI on the effective autoregulation pressure range (Fig. 4) can be similarly rationalized. Under low perfusion pressure, both with and without MVI, the metabolic regulation is fully activated to maximize the vessel diameters and the associated flow (Fig. 5, A and G). Under increasing perfusion pressure without MVI, the vessel diameters rapidly increase under their high transvascular pressure. Thus, the vascular metabolic activation is not needed and is therefore shut down (to around zero level; Fig. 5, K and L), resulting in increased flow and loss of ability to autoregulate. With MVI, even under high perfusion pressure, the vessel transvascular pressure is low, and the metabolic activation is distributed within the limits of 0 to 1 (Fig. 5 D) to retain the ability to autoregulate.
Appendix A: Network flow analysis
Flow in a single vessel
Flow in a single vessel (Fig. 1 A) is governed by the Poiseuille relation as:
where Q(t) is the flow rate, R(t) is the vessel radius, (Pin(t) – Pout(t)) is the longitudinal input/output pressure drop, µ(R) is the dynamic viscosity which is a function of vessel radius, and L is the length of the vessel (Jacobs et al., 2008). The shear stress, τ, is given by
Flow in each vessel is simulated by a validated lumped three-element Windkessel model (Jacobs et al., 2008). Two nonlinear resistors (, ) in series are connected in parallel to a capacitor C (Fig. 1 B). The capacitive element represents the pressure-induced volume change in each elastic vessel. The resistance ℜ (and conductance G), and the capacitance C of each vessel are
The junction of the three elements in the single vessel (Fig. 1 B) is the geometric center of the vessel with an unknown pressure, Pmid. The bifurcation node between vessels is the junction of three resistors, each belonging to a different vessel (Fig. 1 C). The nodal pressure at the junction of the three vessels is Pnode. The unknowns Pmid and Pnode are solved based on conservation of mass—there should be no net volume imbalance at the internal midpoint in the Windkessel element or at the bifurcation node.
Assumptions of the model are (a) each vessel is of uniform cross section and wall thickness, (b) flow in the vessels is laminar, (c) blood viscosity in each vessel is diameter-dependent following measured data of the apparent viscosity in microvascular beds (Pries et al., 1994), (d) active and passive vessel properties are homogeneous in each vessel but vary between vessels, (e) the dynamic extravascular pressure PT is constant along each vessel but varies between vessels depending on their transmural location, and (f) the active vessel response depends on the time-averaged pressure, flow, and metabolic signal.
Network flow analysis
Flow in the multiple vessels network is determined via the iterative solution of a system of ODEs based on the conditions of conservation of mass. It requires that the net flow in each node be zero, and hence for a vessel midpoint,
where is the given input signal of the extravascular pressure, which depends on the myocardial transmural wall location. The mass conservation at the midpoint in each vessel is given by
where Pin(t) is the vessel input pressure signal and Pout(t) is the vessel output pressure (Fig. A1).
The net flow at a bifurcation or trifurcation between a mother and daughter vessels at a designated “network node” is zero. By application of mass conservation at each node, additional equations for the nodal pressures are obtained. For an ith vessel, which is neither source nor sink, mass balance at the vessel inlet (Fig. 1 A) yields
Hence, the inlet nodal pressure as a function of neighboring vessel pressures and conductance is given by
Similarly, applying mass balance at the outlet node of the ith vessel gives
For a source vessel, the inlet pressure Pin(t) and for a sink vessel, the outlet pressure, Pout(t) are prescribed boundary conditions. The governing equations are assembled for the network into a system of ODEs which is written in matrix form as
Expressions for the A and B matrices are given in Appendix B. The trifurcation node is the junction of four resistors, each belonging to a different vessel. Similar to a bifurcation, the nodal pressure (Pnode) at the junction of the four elements is an unknown. The network flow was solved based on mass conservation at each vessel midpoint and at each vessel nodal junction.
The terminal arteriole flow
The flow in the precapillary arterioles, qterm, represents the cardiac perfusion. The flow must satisfy the myocyte demand for oxygen as represented by the target flow qtarget. This is an input to the model that depends on the heart metabolic demand. Several levels of qtarget were used to represent a range of physiological activity levels.
The solution of network flow
Because of the vessels’ elasticity and their interaction with their embedding myocardium (MVI), the flow simulation is highly nonlinear. Network flow solution was thus an iterative solution of the system of ODEs (Eq. A10) subject to the boundary conditions of each case (Fig. A1). The matrices A and B were modified after each iteration. Iterations were continued up to satisfaction of convergence and periodicity conditions to within specified tolerances.
The numerical framework first solves the passive network flow, followed by that of the actively regulated one. Two different schemes were used to solve the actively regulated flow. In cases where metabolic regulation was absent (i.e., only myogenic and/or flow mechanisms active), the solution for a passive vessel network was used as the initial guess to adjust each vessel diameter following the pertinent model equations, according to its predicted pressure and flow rate. When the metabolic regulation was active, the solution from passive network flow was iterated under a genetic algorithm search for the distribution across all the j terminal arterioles of the metabolic signals (Eq. 6), which yield terminal flow rates close to qtarget, up to within a specified tolerance.
The solution of a reference case was performed with parameters listed in Table 1. It served as a baseline for the sensitivity analysis to compare predictions under a range of parameter levels.
The accuracy of numerical flow solution
Simulations were performed to verify the flow periodicity condition, i.e., the smoothness of the transition between nodal pressures from the end of one cardiac cycle to the start of the next. The smoothness tolerance was set to 0.075 mmHg. The convergence of the computational results was estimated based on the network flow solution, where the net inflow/outflow deviation at each time point (from the requisite zero level) was calculated at each vessel midnode and at each coronary bifurcation. This was performed for both the passive and regulated network flow and under both steady as well as dynamic flow conditions. Convergence was satisfied when all flow deviations were <1% of the inflow to the respective node and vessel junctions.
Appendix B: The A and B matrices ()
The network structure matrix in Eq. A10 for a subtree (n = 400; Fig. S1 and Table S1) was used to build the coefficient matrices A (n × n) and B (n × 1). The indices of the nonzero elements at an ith row in A correspond to the ith vessel properties and its neighbors. A vessel connected to bifurcating vessels at its origin and end has indices (i, n-1), (i, n-2), (i, n1) and (i, n2) for the mother, sister, and two daughters, respectively (Table S1). The corresponding conductances are . A vessel with a trifurcating branch at its end has an additional daughter with index (i, n3) with conductance . A vessel with a trifurcating vessel at its origin has an additional sister with index (i, n-21).
For a vessel, i, connected to bifurcations at both its ends, the elements of matrices A are
where Pmid is the pressure in the vessel midpoint, and C is its capacity (Eq. A4).
For a vessel, i, connected to a bifurcation at its origin and a trifurcation at its end, the elements of matrices A are
For a vessel, i, connected to a trifurcation at its origin and a bifurcation at its end, the elements of matrices A are
If an ith vessel is terminating, the elements at each row in A are
For the source vessel, the elements of A are
For all interior vessels, the elements of B are
For a terminal vessel, i, the elements at each row in B are
For the source vessel, the elements of B are
Pin and Pout are respectively the network input and output pressures.
Appendix C: Passive vascular mechanical properties
Vascular tethering tension
Coronary vessels are tethered to the surrounding myocardium by a network of short collagen struts (Borg and Caulfield, 1979; Caulfield and Borg, 1979). These struts prevent the vessels from collapse under negative transvascular pressure that occurs in the deeper myocardial layers (Kajiya et al., 2008). The importance of tethering is clearly seen when considering the wall tension balance under negative transvascular pressure in the absence of tethering. In the case of passive vessels, the tension balance is expressed by
The compliant untethered vessel wall can sustain only a very small level of negative pressure before collapsing. In the presence of active regulation, the wall smooth muscle cells contract, thereby adding to the vessel tendency to collapse. Tethering prevents collapse by adding a positive pressure–like term to the combined tension balance equation, namely,
Two assumptions were adopted in the analysis. First, under negative transvascular pressure (ΔP < 0), the passive wall tension Tpas vanishes. Second, under positive transvascular pressure, the tension in the tethering struts increases quadratically with the gap between the zero-pressure radius (R0) of the myocardial “tunnel” that embeds the vessel and that of the vessel (Rreg). This nonlinearity in the tension gap relationship reflects the gradual recruitment of the nonuniformly undulated struts with stretch. Second, under positive transvascular pressure (ΔP < 0), the tension in the tethering struts was assumed to vanish while the vessel adheres to the myocardial tunnel. The struts tension is thus Cstr (R0– Rreg)2. This expression cannot be used for coronary vessels of all orders with the same value of Cstr because it assumes that the struts density on the myocardial tunnel wall is independent of R0. It is likely that the struts density is higher for smaller vessels with smaller perimeters and vice versa for larger vessels. Hence, the strut pressure–like stress per unit myocardial tunnel area can be expressed as
The strut density on the vessel wall under negative transvascular pressure (i.e., when Rreg < R0) is higher than that on the myocardium by a factor of R0/Rreg, so that the strut stress on the vessel wall is . Finally, the contribution of the strut pressure–like stress to the vessel wall tension (“tethering tension”) is given by
Appendix D: Vessel wall active tension and dynamic stiffness
The radius of the active regulated vessel varies dynamically under a cyclical transvascular pressure. The vessel response is determined by the active dynamic vessel wall stiffness. In vitro studies of isolated vessels (Halpern et al., 1978) have shown that the vessel dynamic stiffness is proportional to the active tension, Tact. This proportionality is in line with findings of in vitro isolated myocytes and myocytes culture (Lipowsky et al., 1978; Campbell et al., 2003; Yadid and Landesberg, 2010). It is thought that in the heart (and skeletal muscle), both the stiffness and active tension are proportional to the number of attached actin-myosin cross-bridges. Because the contractile machinery in VSMC is actin-myosin as well, it is reasonable to adopt the same relationship for the arterial wall, and with the same proportionality constant. To this end, the total vessel wall tension is taken as the sum of the passive and active components. Under equilibrium, this wall tension balances the contributions of the trans-vascular pressure and of the tethering tension. Hence, we obtain
The tension of the passive wall, Tpas, stems from the passive elements. From Eq. 3, it is given by
where Ap, Bp, Cp, and φp are passive vessel parameters (Bp = 0).
The active tension is taken to depend solely on the vessel circumference (or radius) and activation. This is a result of the dependence of VSMC active tension in the vessel wall on the actin/myosin overlap, which is directly related to the wall circumference or radius. Tethering does not affect the active tension–radius relationship. Hence, at a given passive tension, the force balance in Eq. D1 under = 0 determines the active tension as a function of both the Rreg and the activation level A (Eqs. G1 and G4).
The active constitutive relationship allows the determination of the passive vascular dynamic stiffness, kpas, as follows:
When the regulated radius (Rreg) is less than the passive zero-pressure radius (R0), the tethering tension becomes effective. The tethering “stiffness” is defined similar to the passive one (from Eq. D3), as
where the minus sign designates the decrease in vessel radius to increase in tethering stiffness. Because the dynamic stiffness of the active vessel to stretch perturbations was found to be linearly proportional to the active tension (Halpern et al., 1978), the level was evaluated from the data in their Fig. 8 A to be
where C1 = 30.6 mm−1 is the slope of the linear regression and C0 = 4.85 kPa is the intercept. The total vascular dynamic stiffness is the sum of active, passive and tethered stiffness components as
With the given expressions for each stiffness term, the total vascular dynamic stiffness can be evaluated from Eq. D6.
Similar to Eq. D3, the vessel compliance under dynamic loading of around is obtained from the vessel wall dynamic stiffness, k, as follows:
Appendix E: The metabolic flow regulation
Early studies on the vasomotor response in the microcirculation observed that dilation spreads over a much larger area than could be accounted by diffusion (Krogh et al., 1922). More recent studies established the predominant role of the endothelium layer in conducting vasodilatory stimulus (Furchgott and Zawadzki, 1980; Emerson and Segal, 2000; Looft-Wilson et al., 2004) via cell-to-cell coupling (Larson et al., 1983).
To establish the specific pathway by which coronary diameter is regulated has been difficult because of redundancies in control pathways, differences between species, conflicting results in different studies, different (at times opposite) effects at rest versus during exercise (Duncker and Bache, 2008), and at times opposite effects on vessels of different sizes (Gorman and Feigl, 2012). In humans, there are additional uncertainties because of inadequate control of the endothelium state in coronaries of patient subjects (Duncker and Bache, 2008). Notwithstanding these difficulties, it was established that coronary local metabolic control is not primarily a result of adenosine, ATP-dependent K+ channels, NO, prostaglandins, and inhibition of endothelin (Tune et al., 2004; Duncker and Bache, 2008). A number of mechanisms for the initiation and conduction of vasodilation have been proposed (Xia and Duling, 1995; Doyle and Duling, 1997; Rivers, 1997; Hoepfl et al., 2002; Murrant and Sarelius, 2002; Budel et al., 2003; Looft-Wilson et al., 2004; Figueroa et al., 2007; Tallini et al., 2007). A specific pathway that attracts attention in recent years proposes that red blood cells (RBCs), which are responsible for carrying oxygen in the blood, may also act as sensors of oxygen and thereby of metabolic supply/demand imbalance (Ellsworth, 2000). ATP was found to be released from RBCs in response to hypoxia and hypercapnia (Bergfeld and Forrester, 1992; Ellsworth et al., 1995). These conditions occur in the capillaries and venules under high metabolic demand when oxygen supply is lower than demand (Collins et al., 1998; Farias et al., 2005; Gorman and Feigl, 2012). Venules are thus optimally positioned to monitor the metabolic state of the tissue (Jackson, 1987; Segal, 2005).
Based on a number of studies, the adenine nucleotides regulation mechanism was proposed (Gorman et al., 2003, 2010; Farias et al., 2005; Gorman and Feigl, 2012): ATP released by RBCs in the venules under high metabolic demand is broken down to its metabolites, ADP and AMP. All three adenine nucleotides are potent coronary vasodilators (Gorman et al., 2003). They bind to P1 (AMP) and P2 (ATP and ADP) purinergic receptors on the ECs (Gorman et al., 2003; Burnstock, 2007), thereby stimulating endothelial synthesis of NO, which interacts with the smooth muscle cells in the vessel walls to relax and dilate the vessels, thus reducing their resistance to flow (Sprague et al., 1996).
The vasodilatory signal is believed to be conducted (CR) across the capillaries (Tigno et al., 1989; Collins et al., 1998) to the ECs of upstream arterial microvessels, probably via ECs gap junctions (Segal and Duling, 1987, 1989; Collins et al., 1998; Domeier and Segal, 2007; Figueroa et al., 2007). The vasodilatory effect of CR is believed to decay exponentially with distance into the upstream arterioles (Hirst and Neild, 1978; Delashaw and Duling, 1991; Xia and Duling, 1995), as embodied in Eq. 6 and simulated with a decay length constant L0 of 1 mm—in the range of the flow path lengths in the simulated network (Table A1). Additional experimental support for the CR is found in studies in which ATP application inside small arterioles, outside capillaries, and inside venules produced retrograde conducted vasodilatory response (McCullough et al., 1997; Collins et al., 1998; Duza and Sarelius, 2003).
In addition to the sustained and decaying CR, Figueroa and Duling, 2008 found that short stimulation of acetylcholine evoked transient vasodilation that spread along the entire vessel length (up to 2,000 µm) without decaying. Our study relates to the steady-state effect of sustained metabolic demand. The characteristics of that nondecaying signal and its functional consequences for the entire network flow under sustained metabolic demand are not yet clear. Hence, this mechanism was not included in the present study.
Importantly, the CR used in the present study is not restricted to the vasodilatory mechanism of adenine nucleotides but is rather a generic framework that represents a range of possible vasodilatory signaling spreading from the capillaries to upstream arterioles (Segal, 2005). On the other hand, extravascular synthesized metabolites such as muscular released ATP, adenosine, NO, and potassium, which may diffuse radially into the arteriolar walls and relax their smooth muscle cells, were not included in the present study because in addition to the aforementioned earlier observations (Krogh et al., 1922), more recent theoretical analysis showed that solely local responses provide insufficient flow regulation (Lo et al., 2003). In addtion, countercurrent exchange by diffusion of vasoactive arachidonic acid metabolites between paired venules and arterioles (Hammer et al., 2001) was not included here because, as the authors indicate in their article, pairing and close alignment tend to be typical of larger- and intermediate-size venules and arterioles. Yet a major portion of resistance to flow, and therefore of flow regulation, resides in the smaller-size arterioles. These smaller vessels are highly affected by the CR signal because of their proximity to the venules.
Appendix F: Longitudinal distribution of vessel properties
Under ex vivo conditions, the compliant isolated vessels collapse. In vivo, the vessel tethering to the surrounding myocardium supports it from collapse. To comply with these observations, the level of Bp was set to zero, and tethering was incorporated into the model. The isolated vessel passive parameters (Ap, φp, and Cp) were reestimated from the ex vivo experimental data on isolated vessels (Liao and Kuo, 1997) for vessel orders 5–7 under the constraint of Bp = 0. Parameters for vessel orders 0–1 were obtained by curve-fitting Eq. 1 to in situ diameter–pressure data (Kassab et al., 1999), and parameters of order 10 were extracted from a sigmoidal fit of ex vivo data (Hamza et al., 2003).
The passive parameters as a function of vessel radius across all orders were obtained by fitting a quadratic curve through all data points of the parameter, Ap, and linear fits to φp and Cp (Table A2). The vessel radius was transformed to the radius under 80 mmHg pressure to be compatible with the in situ data (Kassab et al., 1993). Given a vessel radius, the fits were used to interpolate the passive material properties between orders 1–8 (Fig. A2).
The longitudinal distribution of myogenic parameters of vessel orders 5–7 was obtained by fitting Eq. 2 to ex vivo data under varying pressures (Liao and Kuo, 1997). Because capillaries (order 0) and large epicardial arteries (order 10) do not exhibit myogenic radius changes, the myogenic amplitude (ρm) of these vessels was taken as zero. The interpolated sigmoidal parameters of other order vessels are listed in Table A3 and graphically shown in Fig. A3. Because of the lack of experimental data for the myogenic amplitude, ρm, for vessel orders 1–5, the myogenic sensitivity curve was extrapolated from order 5 down to order 1 vessels assuming a constant shape factor m = 2, the level estimated from the in vitro data (Young et al., 2012). The myogenic sensitivity, ρm/R, was assumed to be the same for all order 1 vessels.
From Eq. 2, φm corresponds to of highest ΔRm, the comparison of the ex vivo parameter estimates with the data of order 6 vessels showed that the myogenic diameter reduction ΔRm reaches peak levels at a , which is different from the estimates of φm. Hence, φm for order 6 vessels was adjusted to the value of at highest ΔRm. Furthermore, because the () relationship is symmetric around φm (Eq. 2), the estimates of Cm were adjusted to maintain this symmetry. These adjustments have an insignificant effect on the fit of Eq. 2 to the data. Similar to the myogenic properties, the shear properties Fτmax and Kτ from ex vivo data (Liao and Kuo, 1997) were interpolated as a function of vascular cast radii and are shown in Table A4 and Fig. A4.
Appendix G: Wall tension as a function of regulated radius for vessels with no tethering
The regulated radius as a function of mean transvascular pressure is given by
where the passive contribution of the mean transvascular pressure derived from Eq. 3 as
The wall tension resulting from passive elements of the vessel wall as a function of the regulated radius is given by
The wall tension resulting from active elements is the difference between the total wall tension and the wall tension caused by passive elements of the vessel wall and is given by
The wall tension contributions from the passive and active elements in the vessel wall are calculated with Eq. G3 and G4, based on the experimental data of Liao and Kuo (1997). Because data are available at only a few pressure values, sigmoidal models of PDRs of different order vessels were used to calculate the wall tensions.
For the case of active tension for vessel with tethering, the constitutive equation was modified to include additional pressure on wall cause by tethering as given by
The regulated radius, , is an unknown variable and was determined by the iterative solution of the force balance equation as
Appendix H: Time-varying vessel radius
For flow analysis, Eq. D7 allows the calculation of the requisite vessel radius, R(t), along the cardiac cycle. In the embedded and tethered microvessel, the radius variations are likely small enough to retain just the first term in the Taylor series expansion of R(t), in this case
where ΔP(t) is the time-varying transvascular pressure along the cardiac cycle.
Appendix I: Oxygen demand and target terminal flow
If the secondary contributions of dissolved oxygen on the total oxygen content are considered negligible, the oxygen mass balance is specified by
where M is myocardial oxygen consumption for a single terminal arteriole, qterm is the flow in the terminal arterioles, HD is the hematocrit, Sa is arterial oxygen saturation, Sv is the venous oxygen saturation, and c0 is the oxygen-carrying capacity of RBCs. Hence, we obtain the following:
The values of c0,HD, and Sa are directly measurable and assumed to be fixed (independent of M). The venous oxygen saturation Sv is a function of the oxygen consumption M, and their relationship is dependent on oxygen mass balance, ATP release and transport, and the effect of sympathetic inputs on myocardial oxygen consumption. It is thus affected by the combined action of a feedback pathway signal that is determined by the level of plasma ATP in coronary venous blood, and by adrenergic open-loop (feed-forward) signal that increases with exercise (Pradhan et al., 2016). Data for this relation have been provided by Farias et al. (2005) and Gorman et al. (Gorman et al., 2000a,b, 2010). Based on this relationship between M and Sv, the flow in terminal vessels qterm can be directly related to the oxygen consumption M.
This research was supported by the United States–Israel Binational Science Foundation (grant 2009029) and the National Institutes of Health (grant U01HL118738).
The authors declare no competing financial interests.
Author contributions: R. Namani, Y. Lanir, and G.S. Kassab made significant contributions to creation and design of the study, interpretation of data, and writing the manuscript. R. Namani and Y. Lanir developed the algorithms, and R. Namani wrote the code, ran the simulations, and collected the data. R. Namani, Y. Lanir, and G.S. Kassab approve the final copy of the manuscript, agree to be responsible for all aspects of the research in the manuscript, and confirm the role of authorship.
Eduardo Rios served as editor.