A microstructurally motivated framework to study autoregulation in the coronary circulation

M.J.E., H.G., S.B., D.A.B., J.D.T. and C.A.F. conceived the study design. M.J.E. and H.G. developed the model. All authors were involved in data analysis and interpretation. M.J.E., H.G. and C.A.F. prepared figures and drafted the manuscript. All authors reviewed and edited the manuscript. All authors approved the final version of the manuscript. All authors meet authorship criteria, acknowledge accountability for this work and confirm that no qualifying contributors have been omitted.

This study was supported by the NIH R01‐HL158723, U01‐HL135842, and R01‐HL139813.

Figure  1 A shows an overview of the framework. Anatomically, we considered three different symmetrically bifurcating coronary trees at three different myocardial depths: subepicardium (subepi), midwall and subendocardium (subendo), which are subjected to an intramyocardial pressure. Each subtree consists of 12 vessel generations. Vessel wall mechanics are governed by a microstructurally motivated model based on a CMM formulation, treating collagen, SMCs and elastin as load‐bearing constituents. SMC activation modulates vessel diameter, and is controlled by myogenic, metabolic and shear‐dependent autoregulatory stimuli and determined through haemodynamic modelling. Our framework is calibrated using extensive literature data and considers the relative contributions of the three autoregulatory stimuli. Overview of the coronary tree autoregulation framework ( A ), consisting of subtrees at three myocardial depths: subepi, midwall and subendo subjected to spatially varying intramyocardial pressure ( p i m ). Vessel mechanics are governed by a microstructure‐motivated wall model incorporating collagen, SMCs and elastin. SMC activation is controlled by myogenic, metabolic and shear‐dependent autoregulatory stimuli, which are calculated through haemodynamic modelling. The model undergoes a three‐stage calibration procedure ( B ) to establish passive and active properties that reproduce coronary pressure–flow autoregulation. Framework calibration entails a three‐stage procedure (Fig.  1 B ). In Stage 1, passive and myogenic vessel response is calibrated using pressure–diameter myography data (Liao & Kuo, 1997 ). In Stage 2, the Gharahi et al. ( 2023 ) model is used to define homeostatic morphometry and haemodynamics of the three coronary subtrees. An iterative process between Stages 1 and 2 is performed until the vessel mechanics and homeostatic quantities (haemodynamics, morphometry and SMC tone) converge. Lastly, in Stage 3, metabolic and shear‐dependent autoregulatory responses are calibrated to match a given coronary pressure–flow relationship (Dick et al., 2018 ). The result of this calibration procedure is a tree model endowed with morphological, haemodynamic, structural and autoregulatory characteristics of the coronary circulation. The microstructurally motivated wall model, which governs the passive and active behaviour for all vessels, is presented here, as well as the model of SMC activation and autoregulatory stimuli. Vessel wall mechanics are described with a CMM formulation, considering the main load‐bearing constituents: elastin, collagen and SMCs (Humphrey & Rajagopal, 2002 ). A strain energy function w is used to describe the relationship between stress and strain tensors in a non‐linear manner for each constituent. CMMs incorporate the microstructural properties such that wall constituents are constrained to deform together but have distinct mechanical properties and stresses (Humphrey, 2002 ). Each constituent i has a strain energy density ( w i ), which is a function of deformation gradient that maps from constituent i ’s stress‐free configuration to the deformed configuration. The total strain energy density for a vessel is the sum of its constituents: w = ∑ i w i . The function consists of passive contributions for elastin ( w p a s s e ) and collagen ( w p a s s c ) , as well both passive ( w p a s s S M C ) and active ( w a c t S M C ) contributions for SMC. Therefore, we have: (1) w = w p a s s e + w p a s s c + w p a s s S M C + w a c t S M C . A neo‐Hookean model is used for elastin (treated as an isotropic constituent with coefficient c 1 ), and a Holzapfel exponential model is used for the collagen fibres (with coefficients c 2 and c 3 ) and the SMCs’ passive response (with coefficients c 4 and c 5 ). Critically for this work, the active SMC response is modelled as: (2) w a c t S M C = A S m a x M R S M C ρ w a l l λ θ + λ M − λ θ 3 3 λ M − λ 0 2 , where A is an SMC activation level ( 0 < A < 1 ), S m a x is the maximum active stress, M R S M C is the mass per unit area of SMCs, ρ w a l l is the density of arterial tissue, λ θ is the vessel circumferential stretch, and λ 0 and λ M are the zero and maximum active tension stretches, respectively. Other important material parameters in CMMs of vascular tissue are: (1) homeostatic pre‐stretches at which each constituent is deposited ( G θ e , G z e , G h c and G h S M C ), and (2) the fraction of the vessel wall occupied by each constituent (mass fractions ν e , ν c and ν m ). Details on elastin, collagen and SMC constituent orientations within the vessel wall, as well as other mathematical forms of the CMMs are summarized in Appendix A and available in greater detail elsewhere (Baek et al., 2007 ). Using membrane theory (Humphrey, 2002 ), the circumferential wall tension ( T θ θ ) of a thin‐walled vessel is calculated from its strain energy density as: (3) T θ θ = 1 λ z ∂ w ∂ λ θ , where λ z is the axial stretch of the vessel. Furthermore, mechanical equilibrium is given by Laplace's law: (4) T θ θ = p t v D 2 , where p t v is transvascular pressure and D is lumen diameter. Note, p t v is the difference between lumen pressure ( p ) and the intramyocardial pressure ( p i m ): p t v = p − p i m . By combining Equations ( 3 ) and ( 4 ), the vessel diameter can be determined. Autoregulation is incorporated within the CMM formulation by adjusting the SMC activation A in response to autoregulatory stimuli. SMC activation is modulated between fully dilated ( A = 0 ) and maximally constricted ( A = 1 ). Following Carlson et al. ( 2008 ), A is assumed to be sigmoidal and to depend on the total SMC stimulus ( s t o t a l ): (5) A = 1 1 + e − s t o t a l , where s t o t a l depends on myogenic ( s m y o ), metabolic ( s m e t a ) and shear‐dependent ( s τ ) stimuli. Furthermore, SMCs maintain a basal tone under homeostatic conditions ( s h ) (Martinez‐Lemus et al., 2003 ). Thus, s t o t a l can be expressed as: (6) s t o t a l = s m y o + s m e t a + s τ + s h . It is important to point out that under the homeostatic state s h contains contributions from each autoregulatory mechanism, even though at the homeostatic state the three different signals (myogenic, metabolic and shear‐dependent) are zero by construction (see below). Since the relative strength of these control mechanisms depends on vessel size (Arciero et al., 2008 ; Goodwill et al., 2017 ; Liao & Kuo, 1997 ; Stepp et al., 1999 ), normalized diameter‐dependent response curves ( 0 ≤ ϕ ≤ 1 ) for the myogenic ( ϕ m y o ( D ) ), metabolic ( ϕ m e t a ( D ) ) and shear‐dependent ( ϕ τ ( D ) ) control mechanisms are defined below. The myogenic stimulus s m y o is assumed to be linearly dependent on wall tension T θ θ (Carlson et al., 2008 ): (7) s m y o = ϕ m y o a m y o T θ θ − T θ θ h T θ θ h , where a m y o is a myogenic scaling coefficient and T θ θ h is the homeostatic wall tension. From Stage 1 model calibration, we obtain a m y o and ϕ m y o at four representative vessel sizes ( D  = 180, 95, 62 and 37 µm), as discussed in Appendix B , using the data of Liao & Kuo ( 1997 ). To extrapolate ϕ m y o to other vessel diameters, we assume no myogenic response (i.e. ϕ m y o = 0 ) for the precapillary arterioles and for D  >500 µm (Dick et al., 2018 ). ϕ m y o is interpolated across all vessel diameters using a power‐law scheme. As discussed in Appendix B , each of the three subtrees has distinct ϕ m y o distributions, denoted ϕ m y o s u b e p i , ϕ m y o m i d w a l l and ϕ m y o s u b e n d o . Changes in MV O 2 are assumed to be matched by proportional changes in flow to meet the oxygen demand (Arthurs et al., 2016 ). Based on this assumption, the metabolic stimulus s m e t a is linearly dependent on the terminal arteriole flow rate, regardless of the underlying physiological mechanisms involved (Namani et al., 2018 ): (8) s m e t a = ϕ m e t a a m e t a q t e r m − q ^ M V O 2 q ^ M V O 2 , where a m e t a is a metabolic scaling coefficient, q ^ ( MV O 2 ) is the target flow rate determined by a baseline metabolic demand and q t e r m is the terminal arterioles flow rate. The metabolic response curve ϕ m e t a assumes that metabolic signals originate in the capillaries and propagate upstream via endothelial gap junctions (Namani et al., 2018 ): (9) ϕ m e t a = e − L p D L 0 , where L p ( D ) is the path length from the midpoint of a vessel with diameter D to the capillaries, and L 0 is a characteristic decay length, assumed to be L 0  = 0.001 m (Namani et al., 2018 ). An increase in wall shear stress τ induces relaxation of SMCs. The shear‐dependent stimulus ( s τ ) is modelled as: (10) s τ = − ϕ τ a τ τ − τ h τ h , where a τ is a shear‐dependent scaling coefficient and τ h is a homeostatic wall shear stress. The negative sign indicates a vasodilatory response from an increase in wall shear stress. Measurements of fractional dilatation induced by shear stress in swine coronary arteries at four vessel sizes ( D  = 180, 95, 62 and 37 µm) from Liao and Kuo ( 1997 ) were used as a proxy for the shear‐dependent normalized response ϕ τ . Similar to myogenic control, ϕ τ is extrapolated to the entire coronary tree assuming negligible shear‐dependent control (i.e. ϕ τ = 0 ) for the precapillary arterioles and for vessel dimeters ≥730 µm (Namani et al., 2018 ). A power‐law scheme is used to interpolate ϕ τ across all vessel diameters. Note that at the homeostatic state A = A h and s t o t a l = s h . The homeostatic SMC stimulus s h and A h are related by: (11) s h = − ln 1 − A h A h . Two examples of normalized diameter‐dependent response curves are presented in Fig.  2 . Once the myogenic responses are calibrated for vessels within the Liao & Kuo ( 1997 ) dataset (180 >  D  > 37 µm, see shaded region), two different assumptions for myogenic responses outside this range are considered. In Fig.  2 A (Low ϕ m y o ), the myogenic responses are assumed to drop to zero, whereas in Fig.  2 B (High ϕ m y o ), the myogenic responses are assumed to remain constant outside the range. Note that the High ϕ m y o profile is not intended to be physiologically realistic, but is used to assess how different response curves may impact modelling results. Results are presented using the Low ϕ m y o profile, unless stated otherwise. A , Low ϕ m y o and B , High ϕ m y o profiles for myogenic response curves. Metabolic ( ϕ m e t a ) and shear‐dependent ( ϕ τ ) response curves are identical in both examples. The grey box indicates vessel diameters in the Liao and Kuo ( 1997 ) dataset. Intramyocardial pressure ( p i m ) depends on a shortening‐induced pressure ( S I P ) and a cavity‐induced extracellular pressure ( C E P ) (Algranati et al., 2010 ). S I P accounts for myofibril contraction during myocyte activation and is assumed to remain constant across the myocardium. C E P is the pressure transmitted from the left ventricle ( p L V ), and is greater in the endocardium than in the epicardium (Arts & Reneman, 1985 ). The distribution of p i m along each subtree j (subepi, midwall, subendo) is modelled as (Sturgess et al., 2024 ): (12) p j i m D = S I P + C E P = S I P + β j D · p L V , where β j ( D ) is a normalized myocardial depth distribution, with β j = 1 at the endocardial surface and β j = 0 at the epicardium surface. S I P and p L V waveforms are taken from Namani et al. ( 2018 ) (Fig.  3 A ). We assign β j ( D ) distributions that anatomically represent coronary trees penetrating the myocardium from the epicardial surface, and approaching capillary beds located at normalized myocardial depths β c a p , j of 1/6, 1/2 and 5/6 for the subepi, midwall and subendo layers, respectively: (13) β j D = β c a p , j L p D L p t o t a l , where L p ( D ) is the path length between the midpoint of a vessel with diameter D to the terminal vessel and L p t o t a l is the total subtree path length. The pulsatile p i m at the capillary depth β j ( D ) = β c a p , j is shown in Fig.  3 A for each coronary subtree, while the distribution of mean intramyocardial pressure ( p i m ¯ ) across subtree generations is presented in Fig.  3 B . A , pulsatile intramyocardial pressures ( p i m ) at a normalized myocardial depth of β j ( D ) = β c a p , j for each coronary subtree, calculated from left ventricular ( p L V ) and shortening‐induced pressure ( S I P ) via Equation ( 12 ). B , distribution of mean intramyocardial pressure ( p i m ¯ ) across vessel generations for each coronary subtree, based on the normalized myocardial depth distribution from Equation ( 13 ). C , inlet and outlet lumen boundary pressure waveforms from Namani et al. ( 2018 ). Vessel fluid–solid interaction is modelled using Womersley's solution for oscillatory flow in deformable tubes (Womersley, 1957 ). This approach solves for flow and pressure under linearized continuum assumptions, as described in detail elsewhere (Filonova et al., 2020 ; Gharahi et al., 2023 ). The continuum assumption allows pressures and flows to be rapidly computed throughout the coronary tree, circumventing the need to resolve the complex physics (Forouzandehmehr & Shamloo, 2018 ). In this work, Womersley's solution is extended to account for a time‐ and spatially varying p i m . Furthermore, a system of equations is constructed to solve haemodynamics throughout the coronary tree. Given the small diameters of the vessel trees in this study, the interactions between red blood cells and plasma become complex leading to highly non‐linear viscosity ( μ ) behaviour, for example Fåhræus–Lindqvist effects (Fåhræus & Lindqvist, 1931 ). To account for some of this complexity while maintaining a continuum approach for solving flow, in this work we have adopted the in vivo viscosity law from Pries et al. ( 1994 ): (14) C = 0.8 + e − 0.07 D − 1 + 1 1 + 10 − 11 D 12 + 1 1 + 10 − 11 D 12 , (15) h f = 1 − H D C − 1 1 − 0.45 C − 1 , (16) μ 0.45 = 6 e − 0.085 D + 3.2 − 2.44 e − 0.06 D 0.645 , (17) μ = μ 0 1 + h f μ 0.45 − 1 D D − 1.1 2 D D − 1.1 2 , where D is the vessel diameter, ( H D ) is the haematocrit level and μ 0 the blood plasma viscosity taken to be 0.001 Pa.s. Intramyocardial pressure is assumed to be a known periodic quantity, uniform along a single vessel of length L (i.e. ∂ p i m ∂ z = 0 ). In the frequency domain, flows ( q ) and transvascular pressures ( p t v ) are related through an admittance matrix ( Y ) (Mackenzie, 2021 ): (18) q = Y p t v = Y p − p i m , where q = [ q 1 q 2 ] T and p tv = [ p 1 t v p 2 t v ] T , and subscripts 1 and 2 designate the vessel inlet ( z = 0 ) and outlet ( z = L ), respectively. Details on constructing Y are presented in Appendix C . For non‐zero frequencies ( ω ≠ 0 ), Y is: (19) Y = π D 2 4 c ρ b l o o d 1 − g 1 − c o s − ω L c + i s i n − ω L c 2 1 + c o s − ω L c + i s i n − ω L c 2 − 2 c o s − ω L c + i s i n − ω L c − 2 c o s − ω L c + i s i n − ω L c 1 + c o s − ω L c + i s i n − ω L c 2 , where c is the pulse wave velocity, ρ b l o o d is the blood density, i = − 1 and g = 2 J 1 ( ∧ ) / ∧ J 0 ( ∧ ) . J 0 and J 1 are Bessel functions of the first kind and ∧ = i 3 2 α W o m , where α W o m = D 2 ω ρ b l o o d / μ is the Womersley number and μ the blood viscosity. For the steady component of Womersley's solution (i.e. ω = 0 ), Y is (Qureshi et al., 2014 ): (20) Y = π D 4 128 μ L 1 − 1 − 1 1 . Thus, Equations ( 18 )–( 20 ) calculate flow and pressure and the ends of a deformable vessel. At any bifurcation, conservation of flow and continuity of pressure states: (21) p m z = L = p d 1 z = 0 = p d 2 z = 0 , (22) q m z = L = q d 1 z = 0 + q d 2 z = 0 , where m and d 1 , d 2 denote mother and two daughter branches, respectively. By applying Equation ( 18 ) to each vessel and Equations ( 21 ) and ( 22 ) at each bifurcation, we construct a linear system of equations: C x = b , where C is a coefficient matrix, x the vector of unknown flows and pressures, and b is a vector containing the prescribed boundary conditions and p i m waveforms. In this paper, pulsatile haemodynamics are only considered at the homeostatic state, which uses the inlet and outlet pressure boundary conditions taken from Namani et al. ( 2018 ) (Fig.  3 C ). Given a fully calibrated model (i.e. after completing calibration Stages 1–3), we can predict autoregulated states following a certain trigger to the system. The algorithm to predict these autoregulated states is shown in Fig.  4 . Starting with the homeostatic tree morphometry, the admittance Y for the steady component of Womersley's solution (i.e. ω = 0 ) is computed for all vessels (Equation  20 ). Then, accounting for the trigger [e.g. change in perfusion pressure ( p i n )], the steady state haemodynamics are calculated throughout the coronary tree. The updated haemodynamics alter the autoregulatory stimuli ( s m y o , s m e t a , s τ ), and consequently SMC activation A . Next, vessel diameters D throughout the tree are calculated using the updated SMC activation (Equations  3 and  4 ), resulting in a new tree morphometry. The process of updating tree haemodynamics and vessel diameters is iteratively performed until convergence, determined when p , q and D change less than a small number ε between two successive iterations. This algorithm outputs the autoregulated haemodynamics, tree morphometry and autoregulatory stimuli. Iterative procedure to solve for the steady‐state coronary haemodynamics for a given perfusion pressure p in For calculating steady‐state haemodynamics at an autoregulated state, a constant capillary resistance R c a p i l l a r y boundary condition, plus a 20 mmHg terminal coronary venous pressure (Pries et al., 1995 ), is applied to the outlets. The capillary resistance is computed as: R c a p i l l a r y = p h t e r m − 20 m m H g q h t e r m , where p h t e r m and q h t e r m are the terminal arteriole homeostatic pressure and flow, respectively, obtained in Stage 2 model calibration (see Appendix B ). All literature data used for model calibration were obtained from left‐ventricular measurements, predominantly from vessels in the left anterior descending (LAD) perfusion territory. Accordingly, the idealized coronary tree developed in this study is intended to be most representative of an LAD microvascular bed. An outline of the framework calibration is presented in this section and discussed in depth in Appendix B . Briefly, the framework was calibrated in three stages (Fig.  1 B ). In Stage 1, passive and myogenic vessel responses at three myocardial depths were calibrated using experimental myography data across four vessel sizes (Liao & Kuo, 1997 ), while enforcing physiological constraints, including realistic constituent pre‐stretches and a target homeostatic smooth muscle activation of 0.5. In Stage 2, homeostatic coronary tree morphology and haemodynamics were determined using an extended Murray's law optimization approach (Gharahi et al., 2023 ). For each myocardial depth, the subtree is symmetrically bifurcating and has 12 generations, with CMM parameters determined in Stage 1 calibration assigned to individual vessels according to their homeostatic diameter D h (Small Arteries: D h  > 140 µm, Large Arterioles: 140 ≥  D h > 80 µm, Intermediate Arterioles: 80 ≥  D h > 50 µm, Small Arterioles: 50 ≥  D h  µm) and myocardial depth. We used 12‐generation trees because when prescribing literature‐based terminal flow rates of 0.001375, 0.0015 and 0.00165 mm 3 /s for the subepicardial, midwall and subendocardial layers, respectively (see Table  B1 ), the 12‐generation configuration provided the best overall agreement with reported flow–diameter and wall shear stress–diameter distributions (see Fig.  6 ). Note, tree generations are not consistent with the diameter‐defined Strahler order scheme defined by Kassab et al. ( 1993 ). The homeostatic smooth muscle activation ( A h ) was calibrated to achieve target vessel wall thickness while maintaining uniform hoop stress across myocardial layers, consistent with experimental observations of transmural differences in vessel wall thickness (Choy & Kassab, 2009 ). The optimization incorporated experimental bifurcation exponent ( ξ ) measurements to capture size‐dependent heterogeneity in vessel branching patterns throughout the coronary circulation. Lastly, Stage 3 calibration determined the metabolic ( a m e t a ) and shear‐dependent ( a τ ) scaling coefficients using experimental pressure–flow autoregulation data while holding myogenic strength constant. Nelder–Mead optimization minimized the error ( E 3 in Equation  B3 ) between model‐predicted and experimental normalized flow responses across different myocardial depths and perfusion pressures. A sensitivity analysis is performed to assess model parameter sensitivity using the approach of Pradhan et al. ( 2016 ). Three indices, X 1 , X 2 and X 3 , were created to quantify the output sensitivity of calibration Stages 1, 2 and 3, respectively. X 1 and X 3 were defined as the percentage change in the error functions E 1 (defined in Appendix B Equation  B2 ) and E 3 (defined in Appendix B Equation  B3 ), respectively, resulting from a 10% change in individual parameter values: (23) X 1 = 100 E 1 ζ E 1 ζ + 0.1 ζ − E 1 ζ , (24) X 3 = 100 E 3 ζ E 3 ζ + 0.1 ζ − E 3 ζ , where ζ is the parameter value and E 1 ( ζ ) and E 3 ( ζ ) are the associated errors. Because the primary output of Stage 2 calibration is vessel diameters, rather than the value of the cost function used for homeostatic optimization (Appendix D Equation  D6 ), we define the X 2 sensitivity metric as: (25) X 2 = 100 N g e n ∑ i = 1 N g e n D i ζ + 0.1 ζ − D i ζ D i ζ , where N g e n is the total number of generations and D i and is the diameter of generation i . Equation ( 25 ) calculates the average change in homeostatic diameter in response to a 10% change in parameter values. Note, the X 2 metric quantifies the impact of parameter perturbations on the established homeostatic tree, rather than modifying the input parameters during the Stage 2 calibration process itself. We assessed the sensitivity of Stage 2 outputs to material inputs and found the resulting tree morphology to be highly robust, with vessel diameters across all generations exhibiting <0.1% change. This present study is a secondary analysis of previously published data; therefore, additional ethical approval was not required. The primary studies utilized for model development and calibration include Liao and Kuo ( 1997 ), Guo and Kassab ( 2004 ), Kaimovitz et al. ( 2008 ) and Dick et al. ( 2018 ). While Guo and Kassab ( 2004 ) explicitly state compliance with national and local ethical guidelines, the remaining studies do not include a formal ethics statement. However, these studies were supported by NIH funding and were thus required to adhere to the NIH Guide for the Care and Use of Laboratory Animals and undergo oversight by their respective Institutional Animal Care and Use Committees (IACUC).

In this section, we present results from model calibration and demonstrate its application through several examples. Model calibration and homeostatic predictions . First, in the  Framework calibration section , we review the three calibration stages and compare the passive, active, homeostatic and autoregulatory properties of our coronary tree to literature data. In the  Homeostatic coronary flow simulations section , we present homeostatic pulsatile haemodynamic simulations at three myocardial depths and in a small artery and small arteriole. Autoregulatory responses to changes in perfusion pressure . The  Autoregulation in response to changes in perfusion pressure section showcases the autoregulatory responses to changes in coronary perfusion pressure, exploring the relative magnitudes of autoregulatory stimuli. In addition, we compare the transmural autoregulatory behaviour and vessel diameter changes with experimental data not used in model calibration. Sensitivity studies . In the  Assessing the relative contributions of the autoregulatory mechanisms section , relative contributions of the three autoregulatory mechanisms in coronary pressure–flow autoregulation are examined. Finally, in the section  Assessing the impact of vascular dysfunction on autoregulation , we assess how changes in passive and active vessel mechanics affect coronary pressure–flow autoregulation. Our calibrated vessel wall model effectively captures experimentally measured passive and myogenic pressure–diameter relationships across all vessel sizes and myocardial depths (Liao & Kuo, 1997 ) (Fig.  5 A – D ). Calibrated model parameters are given in Tables  1 , 2 , 3 and Fig.  5 E . The passive responses (dashed lines) are independent of myocardial depth, whereas the myogenic responses (continuous lines) exhibit clear depth‐dependent differences at p t v  > 100 mmHg (Fig.  5 A – D ). For each vessel class, the subendo layer has the greatest mass fraction of SMCs (Fig.  5 E ). Notably, large and intermediate arterioles have the largest SMC mass fractions, consistent with having the most pronounced myogenic response (Liao & Kuo, 1997 ). Calibrated passive and myogenic responses across four sizes of coronary vessels [ A : Small Artery ( D h  = 180 µm), B : Large Arteriole ( D h  = 95 µm), C : Intermediate Arteriole ( D h  = 62 µm), D : Small Arteriole ( D h  = 37 µm)], compared to experimental data (Liao & Kuo, 1997 ). The y ‐axis shows diameter normalized by the passive diameter at 100 mmHg ( D h / D 100 ). E , estimated mass fractions of elastin ( ν e ), collagen ( ν c ), and smooth muscle cells ( ν S M C ) vary across vessel sizes and myocardial depths. Constrained mixture model material parameters that are constant across vessel sizes Constrained mixture model autoregulatory parameters that are constant across vessel sizes Constrained mixture model parameters that are variable across vessel sizes * υ m = 1 − υ e − υ c Homeostatic optimization results are summarized in Fig.  6 . Vessel diameters range from 11 to 500 µm across 12 tree generations and show minimal dependence on myocardial depth (Fig.  6 A ). Consequently, diameter bifurcation exponent ξ is also not depth dependent (Fig.  6 B ) and is qualitatively aligned with LAD measurements by Kaimovitz et al. ( 2008 ), though our model predicts larger ξ exponents for D  > 65 µm. The vessel length–diameter relationship matches LAD data presented by Guo and Kassab ( 2004 ) and does not differ between subtrees (Fig.  6 C ). The three tree morphologies have little difference in their diameter vs . generation ( A ), diameter exponent vs . diameter ( B ) and vessel length vs . diameter ( C ) distributions. However, there is a layer dependence on vessel thickness ( D ). The homeostatic SMC activation ( E ) was found to vary between the three layers when targeting equivalent hoop stress distributions ( F ). Mean pressure ( G ), flow ( H ) and shear stress ( I ) vary slightly because of different terminal flows being prescribed for the three tree layers. The grey box indicates vessel diameters in the Liao and Kuo ( 1997 ) dataset. Literature data are presented as mean ± SD. Subepi tree thickness (blue line in Fig.  6 D ) was calibrated using LAD experimental data by Guo and Kassab ( 2004 ) through calibration of the SMC homeostatic activation parameter A h (Fig.  6 E ). This produced a hoop stress distribution, which was then assumed to be identical across all three subtrees (Fig.  6 F ). Using this assumption, distributions of thickness and SMC homeostatic activation A h could be determined for the midwall and subendo trees (green and red lines in Fig.  6 D and E , respectively). Subepi vessels are thicker than subendo vessels, consistent with anatomical measurements in the left ventricle (Choy & Kassab, 2009 ). A h was highest in large vessels and lowest in small vessels (Fig.  6 E ). Of note, the vessel sizes in the Liao and Kuo ( 1997 ) dataset used in Stage 1 calibration span generations 4–7, where noticeable changes in A h occur. Additionally, A h was lowest in the subendo layer, while the subepi and midwall layers exhibited similar activation levels. The hoop stress distributions across all three coronary layers demonstrate fair agreement with Guo and Kassab ( 2004 ) LAD measurements (Fig.  6 F ). Pressure and flow distributions are presented in Fig.  6 G and H , respectively, and are comparable to the left common coronary artery pressure and flow distributions calculated by Kassab et al. ( 1997 ). Most of the vascular resistance is in arterioles with D h  < 100 µm, as shown by the rapid decrease in pressure in Fig.  6 G . Wall shear stress increases roughly 80‐fold across the coronary trees (Fig.  6 I ), with values calculated by our model tending to be greater than that reported by Stepp et al. ( 1999 ) in canine left ventricle epicardial vessels, but in fair agreement with the mathematical model presented by Cornelissen et al. ( 2002 ) based on porcine data. Collectively, Fig.  6 shows that our homeostatic optimization framework produced a coronary tree that is morphologically and haemodynamically consistent with what has previously been reported in the literature. Calibrated values of metabolic a m e t a and shear‐dependent a τ scaling coefficients are presented in Table  2 . Furthermore, the autoregulated, fully dilated ( A = 0 ) and maximally constricted ( A = 1 ) flow–pressure relationships are presented in Fig.  7 A , with the ratio of total flow ( q ) over homeostatic flow ( q h ) in the coronary arteries plotted over varying inlet perfusion pressures p i n . At p i n ∈ ( 80 − 140 ) mmHg, our model fits the experimentally reported flow–pressure relationship well. However, our model underpredicts flow at p i n of 40 and 60 mmHg. Note, in this pressure range the predicted autoregulation curve is closely aligned to the fully dilated pressure–flow relation, indicating a full dilatation of the resistive vessels. To assess how our model responds to changes in metabolic demand, we computed flow–pressure relationships in response to a 50% increase [ q ^ ( MV O 2 ) = 1.5 q h ] and 50% decrease [ q ^ ( MV O 2 ) = 0.5 q h ] in target coronary flow (see Equation  8 ). As expected, the flow–pressure curves show upwards and downwards shifts in response to increases and decreases in metabolic demands, respectively (see Fig.  7 B ). A , the calibrated autoregulated pressure–flow relation for the entire tree (continuous line) compared to experimental measurements (Dick et al., 2018 ), with q h being the homeostatic flow rate. B , a 50% increase or decrease in myocardial metabolic demand shifts the coronary pressure–flow relationship vertically up or down, respectively. Literature data are presented as mean ± SD. Coronary flow and transvascular pressure waveforms are shown for a small artery ( D h  = 246 µm, circular marker in Fig.  8 A ) and small arteriole ( D h  = 25 µm, square marker in Fig.  8 A ) for subepi, midwall and subendo subtrees. Flow and pressure pulsatility is greatest in the subendo layer for both vessel sizes (Fig.  8 B – E ), attributable to the intramyocardial pressure ( p i m ) waveform having the greatest pulse pressure in the subendo layer (Fig.  3 A ). Flow is diastolically dominant for both vessel sizes in the midwall and subendo subtrees. However, blood flow is more evenly distributed throughout the cardiac cycle in the subepi layer, with 53.3% and 47.8% of total flow occurring during systole for the small artery and small arteriole, respectively (Fig.  8 B and C ). Retrograde flow occurs in all subtrees for the small artery, with the subendo layer having the most retrograde flow and the subepi layer having the least (Fig.  8 B ). Conversely, the small arteriole exhibits no retrograde flow in the subepi and midwall layers, and only minimal retrograde flow in the subendo layer (Fig.  8 C ). Notably, the shape of the subepi transvascular pressure waveforms differs from the midwall and subendo waveforms during systole and early diastole (Fig.  8 D and E ), which ultimately impacts flow dynamics. Graphical representation of the three coronary subtrees ( A ). The circular and square markers indicate the location on the subtrees from wherethe pulsatile flow ( B , C ) and transvascular pressures ( D , E ) are taken. Pressure, diameter, wall shear stress, autoregulatory stimuli and SMC activation ( A ) are presented for varying perfusion pressures in a large arteriole ( D h  = 110 µm; Fig.  9 A – G ) and small arteriole ( D h  = 25 µm; Fig.  9 H – N ). In the large arteriole, lumen pressure increases linearly ( A ), whereas vessel diameter ( B ) and wall shear stress ( C ) exhibit non‐monotonic responses to changing coronary perfusion pressure ( p i n ). The myogenic stimulus ( D ) is greater than the metabolic ( E ) and shear‐dependent ( F ) stimuli, indicating myogenic control plays a dominant autoregulatory role in the large arteriole. SMC activation curve ( G ) is steepest in the subendo layer and shallowest in the subepi layer, indicative of a layer‐dependent autoregulatory range. In the small arteriole, lumen pressure increases non‐linearly ( H ), and vessel diameter ( I ) and wall shear stress ( J ) exhibit a response qualitatively similar to the larger arteriole. Contrary to the large arteriole, the small arteriole myogenic stimulus ( K ) is weak compared to the metabolic stimulus ( L ). Shear‐dependent control plays a negligible role in the small arterioles ( M ). The small arteriole SMC activation curve ( N ) exhibits a layer‐dependent response. Lumen pressure ( p ) increases linearly with perfusion pressure ( p i n ) in the large arteriole (Fig.  9 A ) and a slight non‐linear pattern in the small arteriole (Fig.  9 H ). For both vessel sizes, p shows minimal variation between the layers. Diameter ( D ) responses (Fig.  9 B and I ) show both vessel sizes modulating their diameters over the autoregulatory pressure range, with vessel diameters decreasing as p i n increases from approximately 60 to 165 mmHg. Outside this range, diameters increase with increasing pressure, consistent with the mechanical behaviour in vessels not actively modulating tone. Notable differences in the diameter–pressure response between the three layers are apparent, with the subendo layer actively modulating its diameter over the narrowest range of p i n . Wall shear stress ( τ ) increases with perfusion pressure in both vessel sizes until approximately 165 mmHg before declining (Fig.  9 C and J ), which is consistent with the sharp increase in diameter. The change in wall shear stress is layer‐dependent; however, it is qualitatively consistent between the layers. Next, we examine the relative contribution to SMC activation A of the myogenic, metabolic and shear‐dependent stimuli. All stimuli are set to zero at the homeostatic perfusion pressure (here, p i n = 100 mmHg; see Equations  7 , 8 and  10 ). Larger deviations from zero therefore indicate larger contributions from that signal. The large arteriole exhibits predominantly myogenic control (Fig.  9 D ), with the myogenic stimulus lowest in the subendo layer for p i n   100 mmHg. Metabolic and shear‐dependent contributions are comparatively small and show minimal layer dependence (Fig.  9 E and F ). In contrast, the small arteriole shows greater metabolic control (Fig.  9 L ), with myogenic control contributing substantially less than in the large arteriole (Fig.  9 K ). In small arterioles, the metabolic stimulus is lowest in the subendo layer for p i n   100 mmHg, whereas myogenic and shear‐dependent stimuli remain relatively uniform across layers (Fig.  9 K and M ). These differences in stimulus between the vessels are a direct consequence of the assumed normalized response curves in Fig.  2 A . Lastly, Fig.  9 G and N show the relationships between SMC activation A and p i n . The steepest patterns are found in the subendo layer for both vessel sizes, indicating that vessels in this layer reach maximum dilatation ( A  = 0) or constriction ( A  = 1) within a narrower pressure range. The ratio of subendo to subepi flow (ENDO/EPI) at varying perfusion pressures is studied next, using data never used for the calibration of our model. Figure  10 A shows ENDO/EPI data in both anaesthetized (Boatwright et al., 1980 ; Kanatsuka et al., 1990 ) and un‐anaesthetized (Canty, 1988 ) canine experiments in the left ventricle. We run our model with two different assumptions for the normalized myogenic response curve (Low and High ϕ m y o , Fig.  2 ). A homeostatic ENDO/EPI ratio of 1.2 at p i n  = 100 mmHg is enforced by setting the terminal arterioles outflow ratio: q s u b e n d o t e r m / q s u b e p i t e r m  = 0.00165/0.001375, see Table  B1 . From this homeostatic setpoint, the Low ϕ m y o profile produces increased ENDO/EPI ratios for p i n  > 100 mmHg, and decreased ENDO/EPI ratios for p i n   100 mmHg, but not for p i n  < 100 mmHg. A , the calculated ratio of flow going into the subepi layer vs . subendo layer (ENDO/EPI) matches the trends reported in the literature. Both Low and High ϕ m y o models show dilatation (Δ D  > 0) in vessels with a homeostatic diameter ( D h )  160 µm, the Low ϕ m y o model predicts constriction, while the High ϕ m y o model results in Δ D remaining nearly unchanged. The grey box indicates vessel diameters in the dataset of Liao and Kuo ( 1997 ). Literature data are presented as mean ± SD. Next, we explore changes in subepi vessel diameters in response to mild ( p i n = 59 mmHg, Fig.  10 B ) and severe ( p i n = 38 mmHg, Fig.  10 C ) reductions in perfusion pressure. Percentage change in vessel diameter (Δ D ) is compared to experiments in canine left ventricle epicardial vessels (Kanatsuka et al., 1990 ) using the Low and High ϕ m y o profiles (see Fig.  2 ). Using the Low ϕ m y o profile, for both reductions in perfusion pressure, vessels with D h   160 µm constrict (Fig.  10 B and C ). In contrast, the High ϕ m y o response produces similar dilatations for vessels D h   160 µm for both reductions in perfusion pressure (Fig.  10 B and C ). This suggests that differences between Low and High ϕ m y o models are greatest for the larger vessels. In Fig.  11 A , we explore the impact on pressure–flow response of six different combinations of autoregulatory mechanisms (see Equation  6 ). In all six combinations, the basal homeostatic stimulus ( s h ) is maintained. The passive ( A  = 0), fully constricted ( A  = 1) and baseline autoregulated response curves from Fig.  7 A are also presented. The metabolic‐only case (i.e. with s m y o = s τ = 0 ) yields a greater autoregulatory response compared to the myogenic‐ or shear‐dependent‐only cases. The myogenic‐only case maintains some autoregulatory capacity, whereas the shear‐dependent‐only case does not autoregulate. For the myo + met case, the changes in coronary autoregulation are minor compared to the baseline autoregulated case, emphasizing that shear‐dependent control has minimal contributions to pressure–flow autoregulation. These results suggest that the relative contributions of the autoregulatory mechanisms follow the order: metabolic > myogenic > shear‐dependent. Lastly, when replacing the Low ϕ m y o with the High ϕ m y o distribution (Fig.  2 B ), myogenic control substantially strengthens the pressure–flow autoregulatory response. Therefore, these responses are greatly dictated by the assumed normalized autoregulatory strength distributions. Impact of autoregulatory mechanisms on ( A ) the coronary pressure–flow relationship, ( B ) transmural flow distribution (ENDO/EPI) and subepi vascular diameter responses to ( C ) mild and ( D ) severe reductions in perfusion pressure. High ϕ m y o is only plotted in A , as it is included in Fig.  10 . The impact of autoregulatory mechanisms on transmural flow distribution (ENDO/EPI) (Fig.  11 B ) reveals that the myo‐only, meta‐only and myo + meta cases exhibit behaviour qualitatively consistent with the baseline case across perfusion pressures. Conversely, the ENDO/EPI ratio in the shear‐only case exhibits a distinct behaviour, characterized by high sensitivity to minor fluctuations in perfusion pressure around the homeostatic pressure. This suggests that shear‐dependent mechanisms alone are insufficient for regulating transmural flow distribution. Mild ( p i n = 59 mmHg, Fig.  11 C ) and severe ( p i n = 38 mmHg, Fig.  11 D ) reductions in perfusion pressures produce changes in subepi vessel diameters (Δ D ) that are qualitatively similar to each other for all combinations of autoregulatory mechanisms considered. For D h  > 160 µm, the myo‐only and myo + meta cases closely follow the baseline curve, whereas the meta‐only and shear‐only cases experience greater constriction. This is indicative of myogenic control being the most influential in determining the vascular response for vessels D h  > 160 µm. For vessels 50 <  D h  < 160 µm, the myo + meta case experiences the greatest dilatation, due to the absence of the constrictive stimulus provided by the shear‐dependent mechanism. Conversely, the shear‐only case demonstrates the most pronounced constriction, as the reduction in perfusion pressure leads to a decrease in flow (Fig.  11 A ) and a consequent loss of shear‐mediated dilatory stimulus. The myo‐only case best follows the baseline curve in the 50 <  D h  < 160 µm range, also indicative of the myogenic response being the dominant mechanism in this diameter range. However, for D h  < 50 µm, the myo‐only case becomes more constricted than the baseline and converges toward the shear‐only curve. By contrast, both the meta‐only and myo + meta cases converge with the baseline case in this diameter range. This behaviour indicates that the metabolic mechanism is the dominant regulator for the smallest vessels, with myogenic and shear‐dependent mechanisms playing a diminished role. A key benefit of a microstructurally motivated model of coronary autoregulation is its integration of microstructural properties and cellular‐level functions within the vessel wall model, enabling investigation into how structural and cellular changes affect coronary circulation function. We demonstrate the utility of our framework by examining two scenarios of microstructural remodelling: (1) stiffening of collagen fibres and (2) impaired SMC contractility. Specifically, collagen fibres are stiffened by doubling the passive material parameter c 2 , while SMC contractility is impaired by halving S m a x (see Table  1 ). The effects of these two scenarios on flow–pressure autoregulation are shown in Fig.  12 A . For reference, the passive ( A  = 0), fully constricted ( A  = 1) and baseline autoregulated response curves from Fig.  7 A are also presented. Analysis on how changes in vascular mechanical properties impact ( A ) coronary pressure–flow relationships, ( B ) transmural flow distribution (ENDO/EPI), and vascular diameter responses to ( C ) mild and ( D ) severe reductions in perfusion pressure. Stiffening of the collagen fibres results in upward and downward shifts of the fully constricted and dilated curves, respectively. However, coronary autoregulation remains relatively intact, suggesting that collagen stiffening constrains the dilatory capacity of the coronary tree without abolishing its autoregulatory function. In contrast, impaired SMC contractility produces a loss of coronary autoregulation accompanied by an upward shift in the fully constricted curve (Fig.  12 A ). This indicates that reduced SMC contractile capacity severely compromises the ability of the coronary circulation to autoregulate. The impact of these mechanical changes on transmural flow distribution (Fig.  12 B ) shows that collagen stiffening produces a flattening of the ENDO/EPI ratio curve following deviations in p i n when compared to the baseline. SMC dysfunction causes the ENDO/EPI ratio to increase for p i n  > 100 mmHg, but shows little change for p i n  < 100 mmHg. Both collagen stiffening and SMC dysfunction reduce the magnitude of Δ D in response to mild (Fig.  12 C ) and severe (Fig.  12 D ) reductions in perfusion pressure. This demonstrates that both vascular dysfunctions limit the ability of the vessel to modulate diameter, as evidenced by the narrowing between the passive and fully constricted curves in Fig.  12 A . Note, changing the passive and active material properties influences the homeostatic morphology of the tree, as evident by the different homeostatic diameters in Fig.  12 C and D .

The study of physiology and pathophysiology in coronary circulation benefits from accurate mathematical models that can account for both the microstructure and physiological function of arterial networks. In this study, we presented a microstructurally motivated coronary autoregulation model that uses a non‐linear continuum mechanics approach to account for the morphometry and vessel wall composition in three idealized coronary trees. Literature data were used to calibrate and test our model. With some modifications, this model can be applied to morphometry‐based coronary trees instead of idealized symmetric trees. Finally, since our model is based on constrained mixture theory, it could be expanded to also study long‐term growth and remodelling in the coronary circulation in response to hypertension, atherosclerosis, etc.

The coronary circulation has been the subject of many modelling studies over the years. A major focus has been image‐based simulations of epicardial coronary haemodynamics to quantify local flow patterns and assess the functional significance of epicardial stenoses (Kim et al., 2010 ; Lissoni et al., 2025 ; Malvè et al., 2012 ; Taylor et al., 2013 , 2023 ). These studies are typically calibrated to a single haemodynamic state (e.g. resting or hyperaemic conditions) and often account for the microcirculation using lumped‐parameter models. Recent studies have coupled epicardial vessels to synthetic microvascular networks to simultaneously simulate epicardial and microvascular haemodynamics (Kim et al., 2023 ; Menon et al., 2024 ; Papamanolis et al., 2021 ); however, these models do not yet capture the regulatory functions of the coronary microvascular. The coronary circulation manifests a complex dynamic response to changes in perfusion pressure over multiple timescales. These timescales span from autoregulatory variations in active SMC tone, taking place within minutes, to slower changes in vessel microstructure (e.g. stiffening, SMC dysfunction) (Kassab et al., 2002 ; Lu et al., 2011 ; Trask et al., 2012 ), to vascular network angiogenesis or rarefaction over the course of weeks to months (Pries et al., 2015 ). Previous modelling endeavours of coronary autoregulation have focused on capturing key short‐term features using primarily lumped‐parameter approaches (Cornelissen et al., 2002 ; Gharahi et al., 2022 ; Namani et al., 2018 , 2020 ; Xie & Wang, 2015 ). These models lack detailed structural and morphometric information on vascular networks, and are therefore not suitable to predict the long‐term pathophysiological alterations in the coronary circulation noted above. To address this need, this work aimed to produce a unified computational framework integrating vascular responses at the short time scales (e.g. autoregulation) as well as long time scales (e.g. vascular growth and remodelling) while including morphometric, haemodynamics and structural (i.e. vessel wall composition) data (Fig.  1 A ). This design goal resulted in the adoption of a microstructurally motivated modelling framework based on constrained mixture theory, which integrates vessel microstructural properties and cellular‐level functions in a non‐linear continuum mechanics framework. To our knowledge, this is the first implementation of a constrained mixture theory to model coronary autoregulation. The main findings of the model are discussed next. In Stage 1 of model calibration, the vascular wall model was calibrated using passive and myogenic pressure–diameter data from ex vivo myography measurements in swine across four vessel sizes (Liao & Kuo, 1997 ). To apply these data to all myocardial layers, we assumed the passive and myogenic response does not vary transmurally across the myocardium for the transvascular pressures tested experimentally ( p tv ≤ 100 mmHg). To date, only a limited number of studies have investigated whether the myogenic response varies transmurally across the myocardium (Essajee et al., 2026 ; Kuo et al., 1988 ; Sorop et al., 2008 ), and it remains inconclusive whether meaningful differences exist (Dick et al., 2018 ). Model parameters include passive and active material properties, constituent pre‐stretches and mass fractions (Tables  1 , 2 , 3 ). The calibrated model effectively captures the passive and myogenic responses for all vessel sizes, at each myocardial depth (Fig.  5 A – D ). For p tv > 100 mmHg, there are clear differences in the myogenic response between the three subtrees, with the subepi layer experiencing a greater myogenic response compared to the subendo layer. The myogenic response varies between the subtrees partly because vessel wall thickness decreases at greater myocardial depths (Fig.  6 D ), which decreases the total mass of load‐bearing constituents. Future studies probing the myogenic response across a broader range of pressures would elucidate whether myogenic control varies throughout the myocardium. The calibrated SMC mass fraction is greatest in the subendo layer and lowest in the subepi layer for all vessel sizes (Fig.  5 E ). The subendo layer probably has the greatest SMC mass fraction to compensate for its thinner vessels (Fig.  6 D ), enabling it to generate sufficient active tension. Elastin mass fraction remained <10% for all vessel sizes considered, which aligns with these vessels being muscular in nature. To date, mass fractions of coronary artery constituents have only been quantified in large epicardial arteries (Chen & Kassab, 2016 ), while the microcirculation remains understudied. Microscopic analysis of rabbit arterioles showed that SMC content gradually decreases from 100 to 30 µm vessels (Rhodin, 1967 ), consistent with our SMC mass fractions decreasing from intermediate to small arterioles. However, a quantitative assessment of the wall constituents was not performed. Future work quantifying the microstructural composition of the coronary microcirculation will enhance our understanding of coronary physiology and pathophysiology (Hayenga et al., 2012 ). In Stage 2 of model calibration, the homeostatic optimization method proposed by Gharahi et al. ( 2023 ) was used to generate three symmetrically bifurcating coronary trees in subepi, midwall and subendo layers. This approach estimates the homeostatic morphometric (e.g. diameters), structural (e.g. thicknesses) and haemodynamic (e.g. pressure and wall shear stress) characteristics of the vascular trees by minimizing the metabolic and viscous energy dissipation under the constraint of mechanical equilibrium. Furthermore, to incorporate bifurcation diameter exponent ( ξ ) measurements into the homeostatic optimization (Kaimovitz et al., 2008 ), a penalty term was added into the cost function used in Gharahi et al. ( 2023 ). Other approaches for determining a homeostatic state utilize measured or assumed reference geometries and estimate constitutive parameters, including prestretches, to match the prescribed homeostatic geometry (Szafron et al., 2023 ; Wu et al., 2021 ). Both approaches are valid and depend on the information available for calibration. The morphological and homeostatic characteristics of the three subtrees are summarized in Fig.  6 . Swine measurements by Guo and Kassab ( 2004 ) were used as targets for the subepi wall thickness, whereas the midwall and subendo thicknesses were calibrated to match the subepi hoop stress distribution. The decision to match the hoop stress distributions of the three subtrees was based on the experimental findings that vessel wall thickness varies transmurally across the myocardium (Choy & Kassab, 2009 ), suggesting that mechanical stimuli play an important role in determining vessel thickness (Choy & Kassab, 2009 ; Choy et al., 2006 ). Vessel thickness was calibrated by tuning SMC homeostatic activation ( A h ). The predicted homeostatic activation level decreases across generations 1–4 and 7–12, while generations 4–7 deviate from this trend (Fig.  6 E ). The decreasing A h across generations 1–4 and 7–12 is probably a consequence of extrapolating constant material properties to vessels outside the diameter range used for vessel wall model calibration (i.e. all small arteries and small arterioles have the same material parameters). The pressure–diameter data from Liao and Kuo ( 1997 ) encompass vessel diameters of 35–180 µm, whereas our tree spans diameters 11–500 µm. Using variable material parameters might be more appropriate and aid in generating a vessel‐size‐independent distribution of A h (Szafron et al., 2023 ), but we chose not to because our current model already contains many parameters and uses non‐linear constitutive models. Our model would benefit from further characterization of active and passive material properties beyond the experimental range utilized in this study. The discrepancy between our calibrated diameter exponent values for vessels D h  > 65 µm (Fig.  6 B ) and those reported by Kaimovitz et al. ( 2008 ) is due to the two competing terms in the homeostatic cost function (Equation  D6 ), where the first term represents the extended Murray's law and the second term incorporates the experimental measurements. In isolation, the extended Murray's law yields a diameter exponent of ∼3 across the tree (Gharahi et al., 2023 ); however, the inclusion of the Kaimovitz et al. ( 2008 ) dataset pulls the diameter exponent towards the experimental values. Thus, the calibrated homeostatic diameter exponent represents a balance between the energy‐based Murray's law and the experimental data. Nevertheless, our pressure and flow distributions showed reasonable agreement with previous modelling efforts (Fig.  6 G and H ) (Kassab et al., 1997 ), indicative of our homeostatic optimization producing a physiological distribution of vessel diameter (and consequently resistance) across the trees. Note that the inlet pressures for all subtrees were assumed to be 100 mmHg based on the study by Mittal et al. ( 2005 ) that showed that pressure does not significantly drop in vessels of 100–1000 µm. Our wall shear stress tends to be greater than that measured by Stepp et al. ( 1999 ) in canines, but agrees better with previous modelling work by Cornelissen et al. ( 2002 ) which is calibrated using porcine data. Collectively, Stage 2 calibration produced coronary trees at three myocardial depths which have morphological and haemodynamic characteristics which reasonably agree with those reported in the literature (Fig.  6 ). Coronary flow simulations were performed by extending Womersley's solution to account for an external time‐varying intramyocardial pressure (Mackenzie, 2021 ). This approach allows our model to capture important coronary flow physics that are neglected in 0D lumped parameter models, including wave propagations and reflections that influence pressure and flow waveforms (Fig.  8 ). Furthermore, each vessel has a longitudinal spatial variable, removing the spatial ambiguity inherent in lumped parameter models. Flow is diastolically dominant for all vessels in the midwall and subendo layer, whereas flow is more evenly split between systole and diastole in the subepi layer (Fig.  8 ). The subepi layer experiences the lowest p i m , and thus the lowest impedance due to vessel compression. Capturing pulsatile flow waveforms across the heart is important for understanding the spatial heterogeneity of myocardial perfusion and plays an important role in oxygen delivery (Sturgess et al., 2024 ). Collectively, our results show that the magnitude and distribution of p i m are important influences on coronary flow waveforms across the myocardium. Importantly, by extending Womersley's solution to account for time‐varying intramyocardial pressure, we demonstrate that this framework can model haemodynamics across the cardiac cycle, in addition to autoregulatory (minutes) and growth and remodelling (weeks) timescales, spanning key temporal scales of coronary physiology. The coronary autoregulation model is inspired by earlier formulations presented in the literature. We used the SMC activation model proposed by Carlson et al. ( 2008 ), which modulates SMC tone as a sigmoidal function in response to a net autoregulatory stimulus ( s t o t a l ). In addition, we formulated the autoregulatory stimuli in terms of their deviation from a homeostatic value, which was motivated by previous applications of CMMs in the study of long‐term active SMC‐mediated adaptations in cerebral vasospasms (Baek et al., 2007 ). The model resulting from Stage 3 calibration effectively captures coronary pressure–flow autoregulation over a pressure range of 80–140 mmHg (Fig.  7 A ); however, deviations from experimental data are apparent at lower perfusion pressures (40–60 mmHg). Over lower perfusion pressures, the model closely follows the fully passive curve, with SMCs being nearly fully dilated (Fig.  9 ). Since even the fully passive curve falls below the targeted flow levels (Fig.  7 A ), inability to capture the flow–pressure relationship between 40 and 60 mmHg appears related to the passive mechanics of our vessels. As previously mentioned, the passive and myogenic mechanics are determined by ex vivo myography data (Liao & Kuo, 1997 ), which do not capture the impact of myocardial tethering. It has been shown that myocardial tether plays an important role in vessel mechanics (Hamza et al., 2003 ; Hoffman & Spaan, 1990 ; Namani et al., 2018 ), and therefore important for capturing the in vivo pressure–diameter response. Future efforts should include myocardial tethering in the vessel wall model. We explored the response of the model to changes in perfusion pressure. For the two different myogenic response curves considered (Low ϕ m y o and High ϕ m y o ), a reduction in coronary perfusion pressure decreases the ENDO/EPI flow ratio (Fig.  10 A ), consistent with the known susceptibility of the subendo layer to ischaemia during hypotension (Algranati et al., 2011 ; Canty, 1988 ; Kanatsuka et al., 1990 ). Increases in perfusion pressure in the Low ϕ m y o case resulted in increasing ENDO/EPI ratios, as reported by Boatwright et al. ( 1980 ). Thus, the model captures a transmural variability in autoregulatory capacity. Changes in ENDO/EPI are partially a consequence of subendo vessels becoming fully dilated or constricted within a narrower pressure range when compared to subepi vessels, illustrated by the steeper SMC activation curves in subendo vessels (Fig.  9 G and N ). Conversely, increases in perfusion pressure in the High ϕ m y o case resulted in non‐monotonic changes in the ENDO/EPI ratio, demonstrating the impact of the assumed profiles presented in Fig.  2 . We also compared changes in vessel diameters (Δ D ) in response to mild ( p i n  = 59 mmHg) and severe ( p i n  = 38 mmHg) reductions in coronary perfusion pressure with those reported by Kanatsuka et al. ( 1990 ) in canines. In the Low ϕ m y o case, we found that vessels with D h  > 160 µm constricted while D h   160 µm to have relatively unchanged diameter in response to a mild and severe reduction in pressure, contrary to the Low ϕ m y o case. Furthermore, the dilation in vessels D h  < 160 µm was found to be similar between the Low and High ϕ m y o cases, indicative of vessels being fully dilated with both ϕ m y o profiles. For both the mild and severe pressure reductions, our model showed a maximum dilatation at a D h  = 72 µm, and a decreasing dilatation for smaller vessels, a trend that differs from the experimental data reported by Kanatsuka et al. ( 1990 ). However, other models such as that of Cornelissen et al. ( 2002 ) produced patterns consistent with our deceasing Δ D trends for smaller diameters, although their model investigated changes in response to adenosine and l ‐NAME infusion in the coronary microcirculation (Jones et al., 1995 ), rather than changes in response to perfusion pressure. Regardless, it is clear that accurate characterization of this behaviour would require further experimental data on the passive and active mechanics of the smallest coronary arterioles. The Stage 1 sensitivity analysis reveals the sensitivity ranking of passive material parameters to be c 2  >  c 3  >  c 4  >  c 1  >  c 5 , when averaged across the three myocardial layers (Table  1 ). The high sensitivity to collagen parameters ( c 2 and c 3 ) aligns with collagen's highly non‐linear behaviour and its substantial contribution to the vessel wall mass (Fig.  5 E ). The passive SMC parameter c 4 is the third most sensitive, while the SMC non‐linearity parameter c 5 is the least sensitive. Although SMCs contribute minimally to the passive mechanics of large elastic arteries (Lin et al., 2022 ), they are believed to play a more prominent structural role in muscular vessels (Leloup et al., 2015 ). Our Stage 1 calibration identifies a low elastin content, supporting that elastin plays a lesser role in the passive mechanics of the coronary microvasculature. This is reflected by c 1 being the second least sensitive passive parameter. Stage 1 was generally sensitive to changes in both passive and active material parameters. In contrast, Stage 2 calibration was found to be relatively insensitive to perturbations in passive and active material parameters (Tables  1 , 2 , 3 ), with an average diameter change of <5% across the three homeostatic trees in response to 10% parameter perturbations. This difference in sensitivity is probably because Stage 1 is evaluated over a wide range of perfusion pressures, which captures the non‐linear mechanical response of the vessel wall. In contrast, the Stage 2 sensitivity function is evaluated at a single homeostatic perfusion pressure, where the material behaviour is more constrained. Through our Stage 3 sensitivity analysis (Equation  24 ), we identified small arteriole pre‐stretches and S m a x as the most sensitive parameters of the model (Tables  1 , 2 , 3 ). This finding is physiologically reasonable, as small arteriole pre‐stretches directly influence stiffness of vessels that provide the greatest resistance to coronary flow, while S m a x plays a central role in determining the level of active tension. By exploring six different combinations of autoregulatory mechanisms, we were able to investigate their relative contributions to coronary autoregulation (Fig.  11 ). We found that metabolic control was the most dominant in pressure–flow autoregulation, followed by myogenic, and then shear‐dependent (Fig.  11 ). Shear‐dependent control was found to only provide a minor contribution, and removing it only marginally impacted pressure–flow autoregulation, consistent with previous studies by Namani et al. ( 2018 ). The qualitative contributions between the three autoregulatory mechanisms agree with fig. 2A of Carlson et al. ( 2008 ) for arterioles supplying skeletal muscles. We found that switching between the Low and High ϕ m y o profile had a profound impact on the pressure–flow relationship, meaning the relative contribution of the myogenic response is strongly dependent on our normalized myogenic response distribution. Thus, this framework provides a means of investigating how longitudinal autoregulatory mechanisms contribute to global coronary autoregulation, which cannot be assessed experimentally. Additionally, by attenuating or removing individual control pathways, the framework can be used to simulate conditions that alter the responsiveness of coronary autoregulatory mechanisms (e.g. impaired myogenic tone or endothelial dysfunction). Lastly, the utility of our microstructurally motivated framework for pathophysiological studies was demonstrated by examining how collagen stiffening or reduced SMC contractility affects coronary pressure–flow autoregulation. We found that collagen stiffening had minimal impact on pressure–flow autoregulation but significant impact on dilatory capacity (Fig.  12 ). This suggests that while adequate coronary flow may be maintained under resting conditions, the stiffer vessels exhibit insufficient vasodilatation during periods of increased demand. Halving S m a x resulted in complete loss of pressure–flow autoregulation, highlighting the critical role of SMC contractile function in maintaining coronary autoregulatory capacity. These examples demonstrate how the proposed framework can be used to investigate how vascular dysfunction impacts coronary autoregulation. The present study has several limitations. First, the experimental measurements delineating the pressure–diameter relationships of myocardial coronary vessels are scarce. We used four available sets of pressure–diameter data for a wide range of subepi, midwall and subendo vessels from 11 to 500 µm. More pressure–diameter data in the coronary arterioles will enhance the accuracy and predictive capabilities of the model. Second, coronary vascular networks were idealized to bifurcating trees with 12 generations, penetrating the myocardium using a simplified anatomical representation (see section  A model for intramyocardial pressures ). A more realistic reconstruction of the coronary network, potentially by utilizing the diameter‐defined Strahler order scheme (Kassab et al., 1993 ) or microvascular imaging (Sturgess et al., 2025 ), could impact local haemodynamics and autoregulatory responses. Third, the contribution of myocardial tethering is not considered in our model, which plays an important role in coronary vessel mechanics. Fourth, we use simplified models of the autoregulatory mechanisms. The mechanisms involved in the coronary autoregulation, especially the role of metabolic control, have been a largely debated issue in coronary physiology (Kiel et al., 2018 ). In this study we simulated coronary autoregulation using the three main mechanisms and their range of influence as defined in previous modelling studies. Our model, however, could be further refined with improvements in our understanding of the interplay between the different mechanisms of coronary autoregulation. Fifth, the framework uses continuum mechanics assumptions to model haemodynamics and vessel wall mechanics. In the smallest vessels in our coronary tree, discrete cellular and extracellular matrix organization may challenge these assumptions and may ultimately require formulations that explicitly represent structural discreteness (e.g. agent‐based models). However, we demonstrate that the CMM formulation can effectively capture the passive and active pressure–diameter relationship down to the small arterioles (Fig.  5 D ), suggesting that strain–energy functions still provide a useful framework for modelling wall mechanics at this scale. Furthermore, we use an in vivo viscosity law to model haemodynamics and account for non‐Newtonian blood rheology (Pries et al., 1994 ). The continuum formulations used in this study, combined with experimentally determined rheology laws, aim to provide an effective bulk viscosity for microvascular flow computations. However, velocity profiles and near wall shear stress lack the accuracy of formulations that explicitly account for interactions between blood constituents and plasma (Forouzandehmehr & Shamloo, 2018 ). Lastly, while our framework is calibrated using swine data, validation incorporates multiple species (e.g. the canine data in Fig.  10 ). Given the inherent interspecies variability in coronary physiology (Tune et al., 2020 ), the model would benefit from a more uniform approach utilizing calibration and validation datasets from a single species.

The myocardium extracts approximately 70–80% of oxygen from coronary arterial blood under resting conditions. Therefore, shifts in metabolic demand (myocardial oxygen consumption, MV O 2 ) are largely met by changes in coronary blood flow (Feigl, 1983 ). Coronary flow regulation is the result of vessel diameter modulation via local control of the active stress in the vascular smooth muscle cells (SMCs) (Goodwill et al., 2017 ). One aspect of coronary flow regulation is referred to as pressure‐flow autoregulation, in which changes in coronary perfusion pressure are met with adjustments in microvascular resistance to maintain a relatively constant coronary blood flow at a given MV O 2 . Coronary autoregulation relies on SMC reactivity, regulated through myogenic, metabolic and shear‐dependent control mechanisms. Myogenic control is the intrinsic response of SMCs to changes in local wall hoop stress (Kuo et al., 1990 ). Metabolic control is a local feedback mechanism where either an increase in local oxygen consumption or oxygen extraction from the myocardial blood supply triggers a vasodilatory signal (Feigl, 1983 ). Shear‐dependent control is a dilatory mechanism mediated via shear‐induced production of nitric oxide (NO) by the endothelial cells. Studies have shown that the relative strengths of these control mechanisms vary throughout the coronary circulation (Tiefenbacher & Chilian, 1998 ). Myogenic control has been observed to be small in precapillary arterioles and large in 50–150 µm diameter arteries (Liao & Kuo, 1997 ). Shear‐dependent control is most active in the arteries and large arterioles (Kuo et al., 1995 ). Lastly, the signal for the metabolic feedback mechanism is initiated in the capillaries, conducted upstream, and gradually decays for vessels >150 µm in diameter (Kanatsuka et al., 1989 ). In addition to the variable strengths of the autoregulatory mechanisms, the coronary circulation contains several other important heterogeneities. Non‐uniform myocardial‐vessel interactions result in subendocardial vessels experiencing a greater intramyocardial pressure than those in the subepicardium (Zhang et al., 2008 ). Furthermore, vessel morphometry varies transmurally across the myocardium, with vessels of the same diameter reported to be thicker in the subepicardium versus the subendocardium (Choy & Kassab, 2009 ). Collectively, anatomical and functional differences across the myocardium result in a transmural variation in coronary autoregulation (Boatwright et al., 1980 ). The inherent heterogeneous nature of the coronary circulation, combined with the challenges associated with obtaining in vivo measurements, has made our understanding of coronary autoregulation limited and mostly phenomenological. Developing computational models that integrate available experimental data on morphometry, biomechanical and physiological responses can enhance our understanding of coronary autoregulation in health, disease and treatment planning conditions. Coronary autoregulation has been the subject of numerous modelling studies over the last 50–60 years. Virtually all such studies have relied on lumped‐parameter (0D) approaches. Liao and Kuo ( 1997 ) and Cornelissen et al. ( 2002 ) investigated the interaction and balance between autoregulatory mechanisms in the coronary circulation. While these studies incorporated the heterogeneity of the microvascular response, they partitioned the coronary tree into four compartments and ignored the interactions between vessels and myocardium. Pradhan et al. ( 2016 ) developed a data‐driven closed loop model of regulation using in vivo data on coronary flow and oxygen extraction in response to exercise‐induced changes in demand and perfusion pressure. The Pradhan model, which represents parallel control pathways using a block‐diagram approach, does not represent explicit spatial features. Namani et al. ( 2018 ) developed a coronary microcirculation model that integrated dynamic effects of flow with myogenic, shear‐dependent and metabolic feedback control mechanisms in subendocardial and subepicardial vessels. Their model considered simple pressure–diameter rules to define the vessel behaviour. Recently, Gharahi et al. ( 2022 ) and Sturgess et al. ( 2024 ) utilized a non‐linear three‐layer lumped parameter network to investigate autoregulation across the myocardium under varying metabolic demands. While all these studies provide reasonable descriptions of autoregulatory responses, their inherent 0D nature complicates their interpretation. 0D models make it difficult to incorporate microstructurally derived information on structure and function on image‐based networks of vessels. Lastly, given that the short‐term regulation of SMC tone and the long‐term vascular growth and remodelling processes, including pathophysiological responses, are highly intertwined (van den Akker et al., 2010 ), we submit that there is a pressing need to develop microstructurally motivated models of coronary autoregulation. Constrained mixture models (CMMs) have been widely applied to describe the non‐linear mechanical behaviour of arterial tissue, accounting for the contributions of the main load‐bearing constituents (e.g. collagen, elastin and SMCs) (Bellini et al., 2014 ; Ferruzzi et al., 2010 ; Reesink & Spronck, 2019 ). This theory provides a formal means to represent mechanical function of a vessel based on the properties and relative mass fractions of its constituents. Using this theory, long‐term growth and remodelling is represented by changes in the composition of the vessel wall. Previous applications of CMMs in vascular mechanics have focused on large arteries where homeostatic stress has been assumed constant. However, wall shear stress and hoop stress have been found to be size‐dependent through the vasculature (Guo & Kassab, 2004 ; Pries et al., 2005 ; Stepp et al., 1999 ). More recently, CMMs have been applied to study flow‐mediated growth and remodelling (fluid–solid growth) in both large vessels (Baek et al., 2007 ; Figueroa et al., 2009 ) as well as networks of vessels (Gharahi et al., 2023 ; Szafron et al., 2023 ). However, CMMs have yet to be used to study autoregulation in vascular networks. The goal of this study is to develop a microstructurally motivated computational framework of coronary autoregulation, utilizing a recent homeostatic optimization model based on a CMM developed by Gharahi et al. ( 2023 ). The framework provides a means to estimate baseline activation of SMCs and homeostatic stresses, while incorporating experimental data on pressure–diameter, tree morphometry and pulsatile haemodynamics. The utility of the proposed microstructurally motivated framework is demonstrated through different examples of short‐term adaptations to changes in perfusion pressure. A key benefit of this framework is that it provides a method to probe relative contributions of the different autoregulatory mechanisms across different vessel sizes and myocardial depths.

The authors declare that they have no competing interests.

Peer Review History