## Abstract

The mechanisms by which the contracting myocardium exerts extravascular forces (intramyocardial pressure, IMP) on coronary blood vessels and by which it affects the coronary flow remain incompletely understood. Several myocardium-vessel interaction (MVI) mechanisms have been proposed, but none can account for all the major flow features. In the present study, we hypothesized that only a specific combination of MVI mechanisms can account for all observed coronary flow features. Three basic interaction mechanisms (time-varying elasticity, myocardial shortening-induced intracellular pressure, and ventricular cavity-induced extracellular pressure) and their combinations were analyzed based on physical principles (conservation of mass and force equilibrium) in a realistic data-based vascular network. Mechanical properties of both vessel wall and myocardium were coupled through stress analysis to simulate the response of vessels to internal blood pressure and external (myocardial) mechanical loading. Predictions of transmural dynamic vascular pressure, diameter, and flow velocity were determined under each MVI mechanism and compared with reported data. The results show that none of the three basic mechanisms alone can account for the measured data. Only the combined effect of the cavity-induced extracellular pressure and the shortening-induced intramyocyte pressure provides good agreement with the majority of measurements. These findings have important implications for elucidating the physical basis of IMP and for understanding coronary phasic flow and coronary artery and microcirculatory disease.

- coronary flow
- stress and flow analysis
- three-dimensional network geometry
- contraction effects

cardiac contraction has significant impact on the blood flow in the coronary distensible vessels. Myocardial contraction reduces the diameters of intramyocardial vessels, thus increasing their resistances to flow. This “systolic impediment” is the subject of extensive research in relation to *1*) the nature and distribution of the extravascular forces (often termed intramyocardial pressure, IMP) exerted by the myocardium on the coronary vessels and *2*) the manner by which these forces affect the flow in the embedded coronary vessels. The resolution of these issues would allow for a deeper insight into the myocardium-vessel interaction (MVI), which is seminal to understanding coronary pulsatile flow (51) and transmural flow distribution.

Interpretation of the spatial and temporal distribution of IMP from direct in vivo tissue pressure measurements (20, 45) is unreliable due to the cardiac motion, the uncontrolled distortion of the tissue microstructure by the pipette tip (70), and the heterogeneity of this pressure in various microcompartments (57). Downey and Kirk (10) and Spaan et al. (58) assumed IMP to be equal to left ventricular pressure (LVP) at the endocardium and to decrease linearly to zero at the epicardium. Krams et al. (36) found that LVP has only a partial role and proposed (35) the “time-varying elastance” concept to account for the full extent of systolic flow impediment. The elastance concept by itself, however, does not explain why epicardial flows are not inhibited to the same degree as endocardial ones (56). Recent studies have suggested that LVP has a significant effect on time-varying coronary flow as well as myocardial contraction (34). Studies of epicardial lymph pressure (65) and coronary arterial pressure and flow (34) suggest that rather than inducing flow impediment, systolic stiffening shields intramural vessels from the effects of LVP during part of systole (57). Another hypothesis (49) attributes IMP to a combination of LVP-derived interstitial fluid pressure and intramyocyte pressure (48) induced by active contraction.

The mechanistic basis by which the extravascular loading affects the flow in the embedded coronary vessels is not fully understood, as well. Downey and Kirk (10) proposed the “vascular waterfall mechanism” by which increased flow resistance stems from the collapse of intramural vessels. This theory, however, cannot explain the large fluctuations in coronary arterial pressure (58) and the observed flow reversal under certain circumstances (71). Furthermore, such collapse has not yet been observed (23). Spaan et al. (58) introduced the “intramyocardial pump” model, which accounts for the phase lag between arterial and venous flows and for the role of vascular compliance. This concept has since been extended in a number of lumped mathematical models (37, 42, 76) to account for variation in vascular resistance and compliance and for temporal/spatial variation of IMP.

Some local models were developed (66, 67) to study the effects of MVI mechanisms on the local vessel diameter at peak systole and diastole. These models, however, were not integrated into a dynamic network flow analysis and thus cannot provide a complete account of the temporal/spatial coronary flow features. Recently, specific MVI mechanisms were conceptually reviewed (71) but were neither modeled nor tested quantitatively. Hence, the mechanisms underlying MVI and their interplay and effect on flow conditions remain unclear.

In the present study, a three-dimensional morphometry-based model of the coronary network was developed in combination with network flow analysis and local micromechanical investigation based on first principles. Several MVI mechanisms were tested against a wide range of local, regional, and global measurements, both qualitatively and quantitatively. These mechanisms are time-varying elasticity, myocardial shortening-induced intracellular pressure, and ventricular cavity-induced extracellular pressure, and combinations thereof. The relationships of these mechanisms to previously proposed concepts are discussed.

## METHODS

Coronary flow analysis was carried out using each of the above MVI mechanisms. A diagram of the simulation platform used to incorporate all flow determinants is presented in Fig. 1. The scheme consists of a number of submodels and computational processes aimed to reconstruct the network anatomy (Fig. 1*A*), to determine the vessels' reference and loaded diameters (Fig. 1*B*), and to analyze, based on MVI mechanisms (Fig. 1, *C* and *D*), the flow in the entire coronary network (Fig. 1*E*). The network flow boundary conditions [see Online Data Supplement I; supplemental data for this article is available online at the *American Journal of Physiology-Heart and Circulatory Physiology* website], rather than being evaluated from a combined heart-vessel model, were adopted from measured data, this in view of inevitable approximations required in such a complex model and the likewise inevitable need to adjust the model parameters to fit the data.

### Network Anatomy

Reconstruction was based on detailed morphometric data (30–33) and was carried out in two stages: reconstruction of microvascular networks, followed by integration into an intramyocardial coronary network (Fig. 2 and Appendix A).

### Vessel Dynamic Diameter

Vessel diameters are of primary significance to the flow (57). Hence, a major objective of the present study was to develop a mechanical model for the in vivo myocardium-embedded coronary vessel, aimed to determine the in vivo dynamic diameter of each vessel, based on a detailed stress analysis. The scope of this endeavor can be appreciated in light of the large number of vessels involved (arteries, capillaries, veins) and the multiple determinants of the in vivo vessel diameter (Fig. 1*B*). In addition, the vessel stress-free (true) reference diameters are unknown. These were determined by solving the vessel/myocardium stress and strain at five loading configurations for each vessel (Fig. 3). Appendix B outlines the model geometry and constitutive properties, the model equations, the system reference configurations, and the manner by which the loaded dynamic diameter was evaluated throughout the cardiac cycle.

### Flow Simulation

Flow in the coronary network was analyzed using measured data of the system's inputs (Online Data Supplement I): inlet, outlet, and LV pressure waveforms, the corresponding waveforms of cardiac activation, and the sarcomere stretch ratio (SSR). Flow in each vessel was analyzed using a three-element windkessel model consisting of two identical nonlinear resistors and one nonlinear capacitor (Fig. 2, *bottom*). This lumped segment flow model was previously validated (27) against a distributive vascular model (11). At each network bifurcation, mass conservation implies that the sum of discharges Q^{jk} should vanish, i.e., (1a) where P denotes the intravascular pressure in each of the three vessels comprising the *k*th bifurcation, and P is the bifurcation pressure. The vessel hydraulic resistance ℜ is calculated from Poiseuille's law, i.e., (1b) where *L*, *D*, and μ are the vessel length, diameter, and blood apparent viscosity, respectively. The latter was taken to vary with diameter, following Pries et al. (47). Nonlinearity of ℜ(*t*) stems from the vessel elasticity, expressed by the dependence of *D* and *L* on the MVI-dependent extravascular loading, as described below. Mass conservation was implemented into the vessel model, and the resulting equation (Appendix C) was applied in flow analysis of the entire network, subject to given boundary conditions (Online Data Supplement I).

### MVI Mechanisms

Three different interaction mechanisms and two of their combinations were examined.

#### Varying elasticity.

Contractility is taken to affect coronary flow through activation-dependent changes of myocardial stiffness (67). To analyze the effect of varying elasticity (VE) alone, we modeled the myocardium as a hyperelastic solid and ignored the effects of LV cavity pressure (LVP; see Online Data Supplement I) and associated extracellular (interstitial) pressure (Fig. 4*A*). Hence, the boundary conditions for this mechanism were *1*) vanishing extravascular pressure (P_{EV} = 0) and *2*) linear dependence of vessel dynamic diameter on activation (Online Data Supplement I); namely, (2a) where *D*_{p} and *D*_{a} are the vessel diameters under passive and fully active myocardium, respectively (Appendix B). This analysis extends the two-time point (diastole and peak systole) study of Vis et al. (68) to the entire cardiac cycle.

#### Shortening-induced intracellular pressure.

Myocytes for the shortening-induced intracellular pressure (SIP) mechanism were modeled in this study as membrane-contained fluid compartments that surround the vessels (Fig. 4*B*). During shortening, their thickening (lateral expansion of their membranes) is due to an internal pressure elevation. This intramyocyte pressure is transmitted to the vessel due to the impingement by the vessels. This mechanism is based on data and a model (48) showing a linear increase in pressure with contractile shortening. Hence, P_{EV} is proportional to myocyte shortening (expressed by their stretch ratio, SSR) through a scale factor, α, with a baseline value of 0.14 mmHg/%shortening, implying a systolic extravascular pressure of 15 mmHg in the intact tissue. This pressure is lower than the measured intracellular pressure in isolated myocytes (48) but in the range of measured extravascular pressure in intact papillary muscle (21). The baseline value of α was varied through a sensitivity analysis (see below). In this mechanism, the vessel response to pressure is unaffected by activation: (2b)

#### Cavity-induced extracellular (interstitial) pressure.

The cavity-induced extracellular pressure (CEP) mechanism relies on the underlying assumption of both the intramyocardial pump (58) and the vascular waterfall (10) models. P_{EV} is assumed to stem from the LVP alone and is taken to vary linearly (20) with transmural position (expressed as myocardial relative depth, MRD) (Fig. 4*C*). As in the SIP mechanism, activation is not considered in this interaction mode, namely, (2c)

#### CEP + VE.

In this mechanism, the extracellular pressure and varying elasticity were combined. The material law of surrounding myocardium was activation dependent (as in VE), and the extracellular pressure was applied to the myocardium cylinder external surface, i.e., 2d

#### CEP + SIP.

In this mechanism, myocytes were assumed to contract within an LVP-derived pressurized interstitium. Hence, the assigned P_{EV} value equals the algebraic sum of extracellular pressure and the shortening-dependent intracellular pressure such that 2e

A combination of both VE and SIP or of all three basic mechanisms was not considered since the extravascular myocytes can be considered as either a solid (as in the VE mechanism) or a fluid (as with the SIP mechanism), but not as both together.

#### Effects of dynamic axial stretch on flow.

Under all tested mechanisms, the instantaneous diameter and length of each vessel were taken to be affected by the dynamic axial stretch (details in Appendix B).

### Test of the Interaction Mechanisms

Each of the five MVI mechanisms was tested on the basis of its prediction's fit to measured data. This comparison was performed both qualitatively and quantitatively against major measured features of phasic transmural flow and perfusion, velocity, intravascular pressures, and vessel diameters. Transmural perfusion distribution was evaluated from the time-averaged predicted microvascular inlet flow to the network in each specific layer, divided by the weight of the network-perfused tissue. To compare predicted cross-sectional average velocity with data, red blood cells were assumed to be in the center of the lumen (63), and their velocity was considered to be twice that of average velocity (i.e., a parabolic flow profile).

### Effect of Vasomotion

The model was first developed for nonautoregulated vessels and compared with respective data. To compare predictions with in vivo (autoregulated) data, each mechanism was additionally simulated under arteriole diameters reduced by 30% from their calculated dilated values, in line with data comparing nonregulated with regulated vessel diameters (38).

### Sensitivity and Robustness Analysis

Model parameters were taken from reported studies, and no attempt was made to optimize them to better fit the data. In case of uncertainty in reported parameters, sensitivity analysis was performed to examine the robustness of the study conclusions to these parameters. These include the degree of transmission of shortening-induced intramyocyte pressure (48, 72) to neighboring vessels as represented by the proportionality parameter α in *Eqs. 2* and *2e*, which was varied between zero and three times its baseline value. Robustness tests also included the effects of *1*) random structure of the reconstructed microvascular network (tested by using a second stochastic reconstruction), *2*) blood viscosity [tested by replacing the Pries et al. (47) relation with that of Lipowsky et al. (41)], and *3*) in situ axial stretch [tested by varying both tethering stretch values (baseline value is 1.2; Appendix B) between 1 (24) and 1.4 (69) and by neglecting altogether the effect of dynamic stretch (Appendix B) in all vessels]. These tests confirmed that the study conclusions are robust to these perturbations (results not shown).

## RESULTS

Major coronary flow features were compared with the model predictions under each of the tested MVI mechanisms. A summary of the results is listed in Table 1 and elaborated below.

### Transmural Perfusion

In the beating heart, measured subendocardial and subepicardial perfusions were observed to be similar (23). This observation was compared with each of the MVI mechanisms at a heart rate (HR) of 90 beats/min (Fig. 5). Results show that only MVI mechanisms that include the effect of LVP, namely, CEP, CEP + VE, and CEP + SIP, predict a subendocardial-to-subepicardial perfusion ratio (endo/epi ratio) close to unity (0.8, 1.0, and 0.9, respectively). Both VE and SIP mechanisms result in an endo/epi ratio of 1.9, significantly higher than that of the data.

### Total Flow Impediment

The impeding effect of myocardial contraction on total coronary flow is well established (16) and is correctly predicted by all tested MVI mechanisms (Fig. 5). An attenuating effect of contraction was also found under constant LVP (43), suggesting an impeding effect of contraction per se (36). Such an effect of contraction is predicted by VE, SIP, and CEP + SIP (assessed by comparing the predicted total flow under the first 2 mechanisms against flow under no interaction and by comparing flow under CEP + SIP against CEP alone, Fig. 5). Under CEP + VE, however, the total flow is slightly (<1%) augmented, in contrast to the above observation. CEP alone cannot comply with the above observation, since under CEP, cardiac contraction does not affect coronary flow directly but only through changes of LVP.

### Effect of Changes in HR

HR elevation was observed (3) to attenuate both total flow and endo/epi ratio in the canine vasodilated coronaries. An elevation from 60 to 120 beat/min is predicted to change both total flow (by −3, +2, −10, −10, and −13% under VE, SIP, CEP, CEP + VE, and CEP + SIP, respectively) and endo/epi ratio (by 1, <1, −37, −16, and −35%, respectively), as depicted in Fig. 5. Hence flow and endo/epi ratio augmentation under all mechanisms that include CEP is in line with the above observation. These augmentations were found by sensitivity analysis (results not shown) to stem from the lower diastolic time fraction under higher HR values (12).

### Endocardial Sensitivity to Changes in Perfusion Pressure

Under vasodilation, reduction of inlet pressure (P_{A}) was observed to decrease the canine endo/epi ratio (4), whereas pressure elevation resulted in the opposite effect (5). Flow analysis under a reduction of 40 mmHg in P_{A} (from +20 to −20 mmHg of the baseline value) predicts endo/epi ratio changes of +2, +2, −25, −18, and −16% under VE, SIP, CEP, CEP + VE, and CEP + SIP at a HR of 90 beats/min, respectively (Fig. 5). Hence, all mechanisms that include CEP are in line with the above observation.

### Subepicardial Intravascular Pressures

Epicardial arteriolar and venular intravascular pressures were previously found to follow aortic pressure waveform (62). This observation is satisfied by SIP, CEP, and CEP + SIP mechanisms in the subepicardium. In contrast, VE and CEP + VE both predict a subepicardial decrease of arteriolar pressures during midsystole (Fig. 6) in a subepicardial arteriole (30-μm diameter).

### Phasic Diameter

Data on phasic diameter changes were compared with diameter predictions for 60-μm diameter vessels (Fig. 7), since this is the largest diameter of a perfusion vessel and larger (penetrating) vessels are absent in the deeper layers. The simulations predict that CEP, CEP + VE, and CEP + SIP mechanisms result in generally good agreement with the data (Fig. 7). In the VE model, arterial and venous systole-to-diastole diameter changes are −5 and +2%, respectively, in all layers, whereas in the SIP model, they are +1 and −1%, respectively. Both results are much smaller than those of the experimental data.

### Velocity Waveforms

A common observed velocity pattern of subepicardial artery is an early systolic slow down, even flow reversal, followed by a late systolic mild flow and a marked diastolic anterograde flow (8, 10). A similar pattern was also measured in subendocardial arterioles (63). Under VE and SIP, predicted arterial subendocardial velocity remains almost unchanged throughout the cardiac cycle (Fig. 8), in contrast to the data. Under CEP + VE, subendocardial midsystolic velocity is lower than early systolic velocity, also in disagreement with the data. Both CEP and CEP + SIP show good agreement with the data waveform (Fig. 8).

Under mechanisms that include CEP (namely, CEP, CEP + VE, and CEP + SIP), but not under the other mechanisms, the predicted velocity profiles in the largest simulated network artery (450-μm diameter) are mutually similar and in agreement with the profiles measured in the large epicardial coronaries (7, 9) (results not shown; the predictions do not discriminate between the 3 CEP-including MVI mechanisms). It should be cautioned, however, that epicardial vessels are typically larger than 450 μm.

## DISCUSSION

With the goal of contrasting several possible MVI mechanisms against experimental observations of coronary flow, we developed a comprehensive modeling framework based on morphometric data and first mechanical principles. The inputs to the model were taken from measured data. The major conclusion of the study is that neither the effects of LVP (the CEP mechanism) nor those of myocardial contractility (either VE or SIP) alone can account for the transmural and temporal distribution of coronary flow. Only a combination of CEP and SIP agrees well with a considerable body of experimental data.

The results in Table 1 clearly point to the deficiency of mechanisms that exclude CEP. Although CEP predictions generally fit the data well, this mechanism alone cannot account for the observed (43) reduced coronary flow under unchanged LVP levels and increased contractility. CEP + VE can account for this observation only under higher than normal P_{A} values (Fig. 5). Furthermore, under CEP + VE, the predicted subepicardial pressure and subendocardial velocity waveforms differ significantly from the measured one (Figs. 6 and 8, respectively). Most notably, only the combination of CEP + SIP mechanisms is in good agreement with all observations.

The present study extends the two-time point diameter analysis of Vis et al. (66) to a dynamic network flow analysis along the entire cardiac cycle. The results in Fig. 7 demonstrate the importance of these extensions. When considering data of diameter change alone, there is poor discrimination between the different tested MVI mechanisms, especially in light of the high variability of data. In contrast, comparison of flow and pressure data as well leads to a more definitive discrimination.

To model the vascular time-varying elastance (40, 69), Vis et al. (68) considered systolic stiffening, which was found to induce significant flow impediment. Such impediment is also predicted in the present study (Fig. 5, total flow under VE). When combining the effects of LVP and systolic stiffening (CEP + VE), however, the predicted subendocardial extravascular pressure is higher than the intravascular pressure. Hence, myocardial stiffening is shown in this model to actually reduce total flow impediment by moderating the effect of extracellular (interstitial) pressure on vessel narrowing. This increased elastance shielding effect over part of systole was proposed and is compatible with experimental observations of Spaan's group (34, 56).

Interestingly, both VE and CEP + VE models predict midsystolic intravascular pressure attenuation in subepicardial arterioles (Fig. 6), in contrast with available data (62). This may stem from the effect of myocardial activation on vessel pressure-diameter relations (Fig. 3*F*): since predicted transluminal pressures are higher in arteries than in capillaries and veins, the diameter change induced by myocardial stiffening is most significant in arteries. Hence, the longitudinal distribution of resistances is changed (the relative portion of arterial resistance is elevated), causing arterial pressure to decrease during systole.

The results show that the SIP mechanism alone cannot account for the measured systolic flow impediment. To yield higher flow impediment at subendocardium vs. subepicardium (Figs. 5, 7, and 8), a larger subendocardial myocytes shortening must occur, which is not supported by data (50).

The above conclusions of the study can be readily rationalized. Systolic flow impediment is induced by reduction of vessel diameter by external radial forces. The vessel is in contact with, and is therefore loaded by, both myocytes and interstitial fluid. LVP is transmitted to the intramyocardial vessels by the interstitial fluid, since the latter is in contact with the LV cavity. Myocytes impose internal pressure on the vessel, because when they contract, they expand laterally (to preserve their fluid volume), acting like a fluid-filled membrane sac. The myocytes are also in contact with, and are loaded by, the interstitial fluid. Their internal pressure is thus equal to the sum of their shortening-induced pressure as measured in isolated cells (SIP mechanism) and of the interstitial pressure (CEP mechanism), and this combination determines the vessel radial loading.

Myocardial stiffening (VE mechanism), on the other hand, is a material property of the myocardium and as such does not produce radial forces required for flow impediment. Myocyte stiffening may resist systolic vessel expansion, but it does not reduce vessel radius to produce flow impediment. On the contrary, it can shield vessels against the action of external forces such as the LVP, as was found to be the case in the present study and as supported by experimental observations (34, 57). Moreover, in deviance from the original hypothesis of the vascular varying elastance concept (35, 70), there is a basic difference between LV cavity vs. the “vascular cavity.” In the former, myocytes are arranged circumferentially around the cavity. Hence, myocyte contraction produces tangential force that is transformed to radial loading (Laplace law) on the cavity fluid. In vessels, on the other hand, myocytes are aligned normal or oblique to the vessel axis. When they contract, no direct radial force is applied on the vessel cavity except for the intramyocyte fluid pressure (SIP mechanism).

### Other Previously Proposed MVI Mechanisms

Westerhof et al. (71) proposed that geometric changes in myocardium are transmitted to the vessel lumen, thus affecting coronary flow. This mechanism is actually incorporated in all the present MVI mechanisms via the vessel-in-myocardium model (Fig. 3), since myocyte shortening is accompanied by myocyte expansion and vessel deflation. Zinemanas et al. (76) used a single-layer lumped network model and attributed the effect of myocytes stress on the extravascular loading to myocardial volume changes induced by capillary ultra filtration and lymphatic drainage. This model ignores the effects of the coronary network structure and cannot address the important transmural flow differences. Smith et al. (54) attributed the extravascular loading to the myocardial average normal stress components in the radial direction as obtained from a finite element stress analysis during cardiac contraction. The coronary network included only the six largest orders of arteries connected directly to the respective veins. The present model extends that analysis to flow in a realistically interconnected network. Kresh et al. (37) analyzed flow in a lumped model of the coronary vasculature using IMP measurements.

The importance of flow analysis in a realistic network vascular structure (rather than using lumped compartments) is that it allows for determination and coupling of vessels resistance and capacitance based on real properties and basic physical principles and for comparison with data in specific vessels, including microvessels (Figs. 6–8). Attempts to represent portions of the network by lumped models must be validated against a distributive model of a more realistic network such as the ones presented in this study.

### Critique of the Model

Available data that could be used for validation are limited to the subepicardial and subendocardial layers (except for perfusion). In addition, data from different sources are subject to considerable biological variability and at times even to disagreement between studies on different species. For these reasons, and since no parameter adjustment was done, validation was carried out against both quantitative data and a spectrum of generally accepted qualitative patterns. Finally, the concept of shortening-induced intracellular pressure is based on just one study (48) and needs to be confirmed.

The data of Toyota et al. (63) were used for model validation even though it may be biased by the traumatic measurement procedure and possible microvessel occlusion by the used microspheres. The data were used because, to the best of our knowledge, these are the only simultaneous data of subendocardial velocity and diameters. The data reliability is supported by the good agreement of the measured subepicardial velocity profiles with data from larger epicardial coronary arteries (8, 10).

Several approximations and simplifications were used in the model, aimed mainly to reduce its complexity and the associated computational cost: *1*) Vessel flow analysis was based on quasi-static Poiseuillian resistance (*Eq. 1b*), which assumes axisymmetric, steady, laminar, and fully developed flow. Axisymmetry results from the circular cross section used (see below). A posteriori testing of laminarity and quasisteadiness is based on the Reynolds (*Re*) and Womersley (*Wm*) numbers, defined as (3) where ν, ρ, and ω denote blood mean velocity, density, and flow frequency, respectively. The highest predicted values of *Re* and *Wm* (at HR = 120 beats/min and P_{A} = 120/90 mmHg) were 460 and 0.15, respectively, under no MVI and no autoregulation, and 50 and 0.15, respectively, under CEP + SIP and tone regulation, thus justifying the above assumptions.

As commonly used, flow was assumed to be fully developed, which may not always be the case in networks with short internodal distances (75). This assumption may lead to underestimation of vessel resistances (75). The study conclusions are expected to be unaffected by this assumption, since all tested mechanisms were based on this same assumption and are therefore likely to be similarly affected. Furthermore, comparison with the data was based mostly on qualitative trends rather than on exact numerical fit.

*2*) Capillary ultrafiltration was neglected, since it was previously found to have an insignificant effect on the flow in physiological conditions (11).

*3*) The “Gregg effect” (15) was not considered. It suggests decreased contractility under lower perfusion. The rational for this simplification is twofold: it is still controversial (71, 72) if contractility is indeed reduced, and in addition, even if included, it is not expected to change the study conclusions, since all but the predictions of Fig. 5 relate to baseline perfusion and, hence, a constant contractility state. Relating to Fig. 5, decreased contractility could influence the effects on flow of VE and SIP but not of CEP, since only the former two are contractility dependent. Hence, if contractility is reduced, the results under low P_{A} (Fig. 5, *top*) for CEP + VE and CEP + SIP would approach those of CEP alone. Thus the predicted trend regarding the effect of changes in P_{A} (and in HR) on coronary flow and its transmural distribution will remain intact.

*4*) Vein geometry was assumed to be circular, even though veins may have noncircular cross section (32). Unfortunately, the variations of major and minor axes with transmural pressure are not well established. This assumption may be reasonable in diastole, where the major-to-minor axis ratio may be close to unity (13). During systole, however, the vein is more elliptical (13). Again, this assumption is shared by all interaction mechanisms tested and is thus not expected to affect the study conclusions.

*5*) Symmetric transmural networks were used to represent both large arterial and venous trees. This approximation preserves the mean network flow (29) based on average values of the tree morphometric and mechanical properties and is common to all mechanisms tested. Hence, it is not expected to impact the study conclusions.

*6*) Vessel trifurcations and higher order junctions were not included. Such junctions account for only 2% (33) of the arterial tree and 14% (32) of the venous tree. Since the majority of resistance resides in the arterial microcirculation proximal to the venous tree, the model predictions should remain unchanged.

*7*) Energy loss at junctions was neglected following Huo and Kassab (26), who showed this effect on flow to be small. In addition, wave reflections at junctions were not included, based on the previous (74) finding that reflection in small coronary arteries is negligible due to the damping effect of blood viscosity.

*8*) Only limited sensitivity analysis of the effects of model parameters was done. The study conclusions are nevertheless robust for the following reasons. First, parameters were adopted from measurements. In cases where there is significant uncertainty in the values of parameters, robustness of the conclusion was specifically validated for each case (see methods). In addition, most of the model-to-data comparisons were done against qualitative features and are inherently less sensitive to specific parameter levels compared with quantitative comparisons. Moreover, the comparison between mechanisms (e.g., CEP vs. CEP + VE) is a rough sensitivity test, because it assigns two values (zero and baseline) to the associated parameters (of VE in the example). Finally, the results comparing the different MVI mechanisms with the data have been discussed and rationalized based on physical principles.

### Clinical Significance

The present study provides a flow analysis framework based on first principles and compatible with a wide range of coronary data. Hence, it facilitates deeper insight into the physical determinants of clinically used indexes (59) of stenosis severity (e.g., fractional flow reserve) and those of subendocardial hypoperfusion in the presence of upstream stenosis (as emerging from Fig. 5). By doing so, it will be possible to establish the dependence of these indexes on perfusion pressure, diastolic time fraction, and the intramural network structure, as recently proposed (55, 64), based on empirical considerations. In addition, the analysis can provide mechanistic interpretation to data obtained from emerging clinical modalities such as cardiovascular magnetic resonance (44), single-wire pressure-velocity probe (53), and associated wave intensity analysis (9, 60), which have been hitherto difficult to interpret in terms of intramural physiological events (52).

### Summary and Conclusions

Detailed coronary network anatomy was quantitatively integrated with analysis of the vessel-in-myocardium mechanics and a distributive flow analysis to shed light on the mechanism of MVI and to compare five such mechanisms against measured data. The results show that no single MVI mechanism, only the combined effect of LVP-derived interstitial pressure and contraction-induced intramyocyte pressure, fits well to a large body of data. The present findings elucidate the physical basis of IMP and have important implications on understanding and analysis of coronary phasic flow. Furthermore, the study has important clinical significance in shedding light on pathological myocardial perfusion abnormalities and in their diagnosis.

## GRANTS

This research is supported by the U.S.-Israel Bi-national Science Foundation Grant 2003-095 and National Heart, Lung, and Blood Institute Grants HL055554-11 and HL092048.

## DISCLOSURES

No conflicts of interest are declared by the author(s).

## APPENDIX A: NETWORK RECONSTRUCTION

#### Microvascular Network

Stochastic network reconstruction was based on statistical morphometric data (30–33) and was done twice (see Online Data Supplement II), once for flow simulations and once for evaluating the sensitivity to the stochastic network structure. The first network (Fig. 2) features 174 segments: 20 arterioles, 123 capillaries, and 31 venules. Each network is supplied by an order 3 [17-μm diameter, (33)] arteriole (microvascular inlet) and two order −3 (30-μm diameter) venules (microvascular outlet). The volume of myocardium tissue supplied by this network (determined from the network dimensions) is 2.4 × 10^{−6} ml. The network capillary density is roughly 2,800 capillaries/mm^{2}, consistent with measured data (56).

#### Arterial and Venous Trees

To significantly reduce the computational load associated with a full-scale network [consisting of millions of vessel segments (28)] but still retain realistic morphometric features, we placed the reconstructed microvascular networks at four representative transmural locations: at subepicardium (myocardial relative depth, MRD = 0.125), midwall (MRD = 0.325 and MRD = 0.625), and subendocardium (MRD = 0.875). These were interconnected and linked with the major epicardial vessels via intramyocardial arterial and venous tree-like networks, taken to be symmetric and dichotomous (Fig. 2) up to order 8 vessels (28). The latter were reconstructed on basis of morphometric data (32, 33) but assigned identical daughter vessel diameters, lengths, and outlet flow conditions at each bifurcation. The MRD of interconnecting vessels was assigned intermediate values, depending on their transmural location. The reconstructed network has 906 segments, representing the flow in four characteristic myocardial layers.

## APPENDIX B: MECHANICS OF VESSEL-IN-MYOCARDIUM SYSTEM

Following Vis et al. (67), a simplified geometry was considered: each vessel is surrounded by a myocardial tissue of circular cross section (Fig. 3, *A* and *E*). The myocardial outer diameter was chosen to satisfy the measured 1:7 vessel-to-myocardium area ratio (1, 56). Both tissues are considered incompressible and hyperelastic, having a common interface. Each vessel diameter was determined by stress analysis based on the force equilibrium equations of the two concentric cylinders (vessel wall and myocardial tissue) under prescribed loading conditions of dynamic axial stretch λ_{z} (defined as the ratio between loaded and unloaded lengths; see below) and transluminal pressures ΔP. The model equations are presented below, followed by analysis of the vessel and myocardium reference configurations.

#### Model Equations

##### Kinematics.

The axisymmetric mappings between each pair of configurations *i* and *i+1* (Fig. 3) in cylindrical coordinates is (*r*, θ, *z*)_{i} → (*r*, θ, *z*)_{i+1}, prescribed by (B1) where OA denotes the opening angle, *L* is the cylinder length, and the stretch ratio Λ_{i+1,i} = *L*_{i+1}/*L*_{i}. The incompressibility constraint implies that (B2) The Green-Lagrange strain is **E** = (**F**^{T}**F** − **I**)/2, where **I** is unit matrix and **F** is the deformation gradient for each mapping between configurations. Assuming no twist, **F** is given by (B3) Explicit expressions for the deformation gradient of each mapping are given in Online Data Supplement III.

##### Equilibrium equations.

The Cauchy stress tensor **T** is derived from the strain energy function *W* of each material (vessel wall and myocardium; see below) via the hyperelastic relationship **T** = −P**I** + **F**(∂*W*/∂**E**)**F**^{T}. The equilibrium equations in the circumferential and axial directions imply that the Cauchy stress components *T*_{rθ} and *T*_{rz} vanish. The radial force equilibrium equation is ∂*T*_{rr}/∂*r* + (*T*_{rr} − *T*_{θθ})/*r* = 0. By applying the axial and radial equilibrium equations, the external axial force F_{z} and the transluminal pressure ΔP can be expressed in terms of the components *T*_{ij} of the tissue stress tensor **T** as follows (25): (B4)

#### Vessel and Myocardium Constitutive Properties

Vessel wall mechanics was studied in the left anterior descending coronary artery under positive transluminal pressures (69) and is represented by Fung-type exponential material law (see Online Data Supplement IV). The corresponding parameters were estimated for 10 samples. The passive and active myocardial mechanics were previously studied on seven samples and described by invariant-based material laws (40).

To obtain a representative vessel-myocardium pair, each of the 10 vessel parameter sets (69) was combined with each of the 7 parameter sets of the passive myocardial samples (40), and properties of each of the vessel-myocardium pairs was used as inputs to evaluate the passive in vivo pressure-diameter relationship. The pair having the best fit to the in situ data of swine large coronary arteries (19) (Fig. 3*F*) was selected for further analysis of MVI mechanisms.

#### Reference Configuration

Stress analysis requires knowledge of the reference stress-free configuration of both vessel and myocardium. These references have to be estimated to determine the true states of stress and strain. The available data consist of statistics of the in situ diameter, length, and wall thickness (17, 30, 32, 33) taken under vasodilation, no myocardial activation, fixed stretch, and intravascular pressure (vessel cast configuration). The data, however, do not account for transmural morphometric heterogeneity. To incorporate the latter, the reconstructed diameters were first modified by up to ±10% from the vessel cast values in a linear transmural manner to comply with the observed twice higher endocardial than epicardial flows in diastolic vasodilated hearts (14). To obtain the reference configurations, we solved the vessel-myocardium equilibrium equations (*Eqs. B4*) subject to the relevant loading boundary conditions. The latter are *1*) the cast configuration (Fig. 3*A*), *2*) the unloaded configuration (Fig. 3*B*), *3*) the untethered configuration (Fig. 3*C*), and *4*) the stress-free (reference) configuration (Fig. 3*D*).

##### Cast configuration.

Kassab and coworkers (17, 30, 32, 33) measured diameters under fixed cast pressure. The loading conditions in this configuration are ΔP^{v} + ΔP^{m} P_{cast}, where the superscripts v and m denote vessel and myocardium, respectively. Each vessel's specific cast pressure (P_{cast}) was taken from a steady-state coronary flow analysis (46) also based on Kassab's data.

##### Unloaded configuration.

The transition to this configuration is prescribed by the mapping from the cast configuration (Fig. 3*A*). The loading conditions in this configuration are F + F = 0; ΔP^{v} + ΔP^{m} = 0. Based on previous data (18), the axial stretch λ_{z} was taken to remain constant during this mapping.

##### Untethered configuration.

Unloaded coronary vessels are not stress free. When myocardial tethering is removed, large epicardial arteries were found to shorten by 0% in human ex vivo (24) and 40% in swine (69). The untethered configuration was obtained by mapping from the tethered unloaded configuration (Fig. 3*B*). The loading conditions in this untethered state are F = 0; F = 0; ΔP^{v} + ΔP^{m} = 0, where it is assumed that the two cylinders maintain a common but stress-free interface.

##### Stress-free (reference) configuration.

The tethered-free vessel and myocardium are still not stress free, but rather loaded by internal residual stress. The magnitudes of these stresses are quantified by the measured opening angles (OA) of the corresponding cylinders when cut open. For coronary vessels, OA^{v} are specimen dependent. For the myocardium, OA^{m} = 2.75 rad (39). The stress-free reference configuration is obtained by mapping between the configurations represented in Fig. 3, *A* and *D*.

Stress and strain analysis is possible if the radii and stretch ratios at each configuration are known. Since they are not, they were solved by applying the vessel-in-myocardium model equations under the six loading boundary conditions listed above (1 in the cast configuration, 2 in the unloaded configuration, and 3 in the untethered configuration), subject to the assumptions of incompressibility of and common interface between cylinders (see Online Data Supplement V). The solution of this highly nonlinear system was obtained by using the MATLAB ga code for genetic algorithm searches and MATLAB fsolve.

#### Loaded Vessel Diameter

With the stress-free configuration determined, the vasodilated vessel diameters (Fig. 3*E*) were evaluated (using the above MATLAB codes) to optimally satisfy force equilibrium (*Eqs. B4*) subject to the vessel internal and external (MVI dependent) loading conditions. This was done for arteries, capillaries, and veins under prescribed conditions of transluminal pressure, dynamic axial stretch, and myocardial activation. To obtain the corresponding in vivo diameters required for the network flow analysis, we needed modifications to account for the vessel dynamic axial stretch during the cardiac cycle and for the autoregulatory myogenic response and myocardial state of activation, as discussed below.

##### Dynamic vessel stretch.

Vessels are dynamically stretched during the cardiac cycle, together with the surrounding myocardium. This stretch λ_{z}, when added to the constant in situ tethering stretch, affects the vessel diameter. This effect can be evaluated by solving the vessel-in-myocardium model at each time point and for each vessel. To reduce computational load, the stretched diameter was evaluated by linear interpolation between the vessel highest and lowest levels of stretch. The maximum and mean errors associated with this interpolation were examined and found to be 4 and <1%, respectively, in all vessels at all pressures, in either passive or fully active myocardium. The above stretch depends on the vessel orientation. Vessels oriented along myocytes stretch in proportion to the myocytes stretch, whereas vessels perpendicular to myocytes stretch in proportion to myocytes thickening and therefore shorten during myocyte elongation. Hence, penetrating vessels (orders 5–8) were assumed to stretch in proportion to the myocardial wall thickening (see Online Data Supplement I), whereas capillaries (except for cross-connections) were assumed to stretch in proportion to SSR (Online Data Supplement I). Other vessels were assumed to be randomly oriented and therefore unaffected by myocardium contraction (i.e., no dynamic stretch).

##### Myocardial activation.

Myocardial activation was accounted for by calculating the loaded diameters under passive and fully active myocardium (*D*_{p} and *D*_{a}, respectively). The diameter under intermediate activation was determined as a linear interpolation between these two states.

##### Autoregulation.

The entire complexity of tone regulation (8) was not explicitly accounted for, yet to compare predictions with in vivo (autoregulated) data, arterioles diameters were reduced by 30% from their calculated dilated values, in line with data for coronary arteries of <200 μm in diameter (38). Hence, comparison to nonregulated data was preferred whenever possible.

#### Computational Scheme

*Equations B1*–*B4* were used to evaluate the dependence of vessel diameter *D* on the cast value *D*_{cast} and on the transluminal pressure ΔP. These values were presented as surfaces of *D* = *D*(*D*_{cast}, ΔP) for 12 cases: one for each combination of vessel type (artery, capillary, vein), lowest and highest axial stretch, and passive and active myocardium. Based on the experimental results of Hamza et al. (19), the above response surfaces were fitted for each vessel by a sigmoid function expressed by (B5) where *D*_{+} and *D*_{−} denote maximum inflation and deflation diameters, respectively, and ΔP_{1/2} is the transluminal pressure corresponding to the average of *D*_{+} and *D*_{−}. These parameters were estimated to best fit the diameter vs. pressure curve of each vessel. The maximum and mean errors associated with this fit relative to the exact solution were examined and found to be 6 and <1%, respectively, in all vessels under all simulated loading conditions.

## APPENDIX C: NETWORK FLOW ANALYSIS

Conservation of mass requires that the difference between in and out discharges of each vessel *n* (Fig. 2, *bottom*), Q and Q, respectively, should equal the time derivative of the vessel's volume (V^{n}), i.e., (C1) where P and P denote vessel inlet and outlet pressures, respectively. The hydraulic capacitance, *C*(*t*), was defined as (C2)

*Equation C2* extends the common definition of capacity in nonlinear intramyocardial pump models (e.g., Refs. 6, 27), *C*(*t*) ≡ dV/d[P_{IV}(*t*) − P_{EV}(*t*)], to the case in which changes in vessel lumen volume are caused not only by transluminal pressure but also by changes in axial stretch and in myocardial activation. Since each bifurcation pressure P_{Bif} equals either P_{in} or P_{out} of the vessels forming that bifurcation, *Eq. 1* can be combined with *Eq. C1*, resulting in a system of *N* nonlinear ordinary differential equations (*N* denotes the number of network vessels), which was iteratively solved using the MATLAB ode15s solver until it satisfying the periodicity condition (see details in Online Data Supplement VI). The solution process requires boundary conditions of network inlet and outlet pressure (Online Data Supplement I) as well as MVI-dependent extravascular loading (*Eqs. 2*).

- Copyright © 2010 the American Physiological Society