Impact of Coronary Bifurcation Morphology on Wave Propagation

The branching pattern of the coronary vasculature is a key determinant of its function and plays a crucial role in shaping the pressure and velocity waveforms measured for clinical diagnosis. However, although multiple scaling laws have been proposed to characterize the branching pattern, the implications they have on wave propagation remains unassessed to date. To bridge this gap, we have developed a new theoretical framework by combining the mathematical formulation of scaling laws with the wave propagation theory in the pulsatile flow regime. This framework was then validated in multiple species using high-resolution Cryomicrotome images of porcine, canine and human coronary networks. Results demonstrate that the forward well-matchedness (no reflection for pressure/flow waves travelling from the coronary stem towards the microcirculation) is a salient feature in the coronary vasculature, and this result remains robust under many scenarios of the underlying pulse wave speed distribution assumed in the network. This result also implies a significant damping of the backward traveling waves especially for smaller vessels (radius < 0.3 mm). Furthermore, the theoretical prediction of increasing area ratios (ratio between the area of the mother and daughters vessels) in more symmetric bifurcations found in the distal circulation was confirmed by the experimental measurements. No differences were observed by clustering the vessel segments in terms of transmurality (from epicardium to endocardium) or perfusion territories (left anterior descending, left circumflex and right coronary artery). The validation of the proposed theory in multiple species reveals the equivalence between scaling laws and well-matchedness in the vasculature. Moreover, it captures the role of pulsatility in optimal vascular designs. This demonstrates the forward well-matchedness of the coronary bifurcations, whereas backward waves are damped asymmetrically at junctions.


Introduction
The branching structure of the coronary network is an important determinant of its function.
Understanding the distribution of flow, volume, resistance and the resulting pressure and velocity waveforms cannot be achieved without accounting for the underlying design of the vascular network. Since the pioneering work of Murray (40,41), the morphometric relationships between branching vascular segments have been codified into different forms of mathematical relationships, or scaling laws (5,13,23,24,26,76,77,78). Whereas Murray's original work considered the vascular network as an energy-minimizing structure balancing flow delivery and metabolic costs, subsequent work has explored uniform wall shear stress (24,76), myocardial mass (5), vascular volume and cumulative lengths (23,76) as determinants underlying the design principles.
The current leading theory in vascular network design principles is arguably the volume-scaling law (HK law) (23,24), which, based on a set of biophysical assumptions coupled cumulative vessel volume and length with the vessel diameter. This approach produced equivalent results to an earlier minimum-energy model (76) that generalized Murray's formulation from a single bifurcation to the whole arterial network. The discrepancies that were found among different vascular beds were explained on the basis of their varying metabolic-to-viscous power dissipation ratio (26).
Quantitatively, these scaling laws serve a useful role in interpreting angiographic data in the context of growth (4), disease and inter-species comparison (26), in uncovering implied relationships between parameters (56) and estimating the distal flow resistance in myocardial blood flow, e.g. for CFD assessment of FFR (64).
On the other hand, the branching parameters measured in real vascular networks exhibit a large scatter (46,77). Moreover, the fitted power law exponents reported in the literature vary significantly depending on the range of vessel diameters analyzed (67). This heterogeneity had been characterized as possessing features of multifractality using human coronary data, leading to a speculation that such a property may endow the tree with an enhanced ability to distribute impedances along its structure (78). The existing scaling laws provide no explanation for the origin of such heterogeneity. Furthermore, there is evidence that the theoretical framework of scaling laws leads to a shallow optimum -that is to say, the cost of departure from the power law relationship is minimal over a very broad range of exponents (57). Combined with high heterogeneity, it undermines the significance of the putative design principles proposed.
Aside from the physical organization of the vessel segments, pulsatility in flow is also a central feature that is strongly characteristic in the coronary circulation in both large (36,55) and small vessels (66), in that the magnitude of the pulsatile component is not necessarily dominated by the steady component. Surprisingly, to date investigations of coronary network scaling laws have been conducted largely independent of the wave phenomena. This is surprising given that it is the generation, propagation and the reflection of the pulse waves that shape the pressure and velocity dynamics throughout the network. In particular the rhythmic ejection of blood by the left ventricle and the cyclic myocardial contractions that squeeze the embedded vasculature are responsible for the proximal and distal generation of pressure and flow waves. These generated waves then propagate forward (from the coronary root towards the microcirculation) and backward throughout the network, where at each bifurcation encountered, partial wave reflections may alter the waveform shapes. The amplitude of the pulse wave measured at any point in the network therefore carries an imprint of the interaction the propagating wave has undergone with the underlying network structure.
With newer research tools such as wave intensity analysis (WIA) gaining momentum in the clinical arena (8,62,63), the clinical diagnostic potential of altered arterial pulse waves in disease has received much attention in recent investigations (7)(8)(9)55). Specifically, the cumulative energy (area of the waves) and peak energy (peak of the waves) have been used in the clinical settings to enhance our understanding of various disease processes affecting the coronary vasculature (e.g. aortic stenosis (7), LV hypertrophy (8)) and to define the impact of various interventions (e.g. biventricular pacing (34), intra-aortic balloon therapy (10)). Most recently, the demonstrated prognostic benefit of WIA in the setting of acute myocardial infarction has enhanced prediction of long-term myocardial recovery (9). However, while it is known that the LV ejection/suction are the principal factor generating the forward travelling waves, the main mechanisms underlying the backward travelling waves are yet to be elucidated, and confounds the conclusion derived by applying cWIA. Specifically, it is proposed that the myocardial compression/expansion on the intramural vessels and the diastolic reduction of resistance at the microcirculatory level are the two main contributors underlying the backward wave origin (8,55). However, which of the two is the dominant mechanism and how this dominance varies depending on the transmural location and the vessel size is still under debate.
At present little is known about how the branching structure of coronary vasculature is linked with the observed wave propagation behavior. A rare, and perhaps the earliest, work on this subject (1) was the first to propose that the branching structure and mechanical properties of the coronary vessels larger than 0.5mm of diameter are well-matched, so as to facilitate the forward-travelling waves to traverse a junction without being dampened significantly. The presented experimental evidences compared well with their theoretically optimal branching structure. However, this theory was based on an assumption of uniform distensibility in branching vessels, which has not been confirmed for consecutive generations of a given vascular network. This work, therefore, as with many other investigations in arterial wave propagation, suffered from difficulties in determining pulse wave speed and thus relied on extrapolation of results measured in larger segments.
Accurate measurement of coronary pulse wave speeds is a challenge that persists even today, and there is a large scatter in the values employed in the literature ranging from 5-10 m/s in earlier studies (1,77), to higher values 15-25 m/s estimated by more recent ComboWire measurements (48,50). In the analysis below, we show that evoking the theory of scaling law can help to overcome this experimental difficulty.

4
In short, the existing scaling laws capture well the broad behavior of the vascular branching pattern in an averaged sense. Heterogeneity, however, is better addressed by assessing the wave reflections (and the encompassing theory) since they are inherently a local phenomenon. As demonstrated herein, the application of scaling laws to the interventionally-relevant epicardial segments results in a large spread that cannot be explained by the existing theory. By unifying the theoretical framework of the scaling law with those of pulsatility and wave propagation, we show that phasic aspects of coronary flow -which, after all is one of its defining characteristics -is an indispensable determinant of its structural design that can contribute to the observed heterogeneity. Accordingly what we propose here is not a new scaling law, but a generalized framework by which to elucidate vascular branching structures.
The rest of the paper is organized as follows -we begin by recapitulating the theories of scaling laws and wave reflection in vascular networks, before developing the new unified model of the organization of branching structure under pulsatile flow. Following this, the physiological range of Womersley numbers is considered to incorporate the frequency-domain behavior into the analysis. Then, experimental validation comprising high-resolution imaging and segmentation of coronary networks from several species are described and applied to the proposed theory. Finally, model results are applied to the new and existing experimental data to explain the spread, to re-visit the well-matchedness hypothesis and to investigate regional and scale-specific differences in the coronary branching patterns. In addition, existing scaling laws that were constructed in the absence of pulsatile flow are reassessed to show that each proposed law implies a specific branching pattern, which has not been evaluated to date in terms of resulting wave propagation.

Methods
We begin with a brief presentation of the mathematical background of the scaling laws regarding the coronary branching pattern (23,24,40,41,78). Subsequently, the one-dimensional blood flow theory is introduced along with the concepts underlying wave reflection at bifurcations. These two theoretical frameworks are then combined together incorporating the impact of flow pulsatility (using Womersley's approach) into the analysis. Finally the experimental protocol and imaging processing pipeline employed for coronary anatomical reconstruction are described, followed by the model validation.

Theory of scaling law at bifurcations
As described in the Introduction, a scaling law is an analytic formulation describing morphometric relationships, often parameterized by the vessel lengths, diameters and arterial volume between a feeding segment and the perfused subtree (13,24,78). Since the impact of a scaling law on wave reflection at bifurcations is of primary interest in this study, we focus on a specific aspect of the scaling laws: the relationship between the mother vessel diameter (dm) and the diameter of the daughter vessels (dd1, dd2) at each bifurcation.
This relationship can be written in a generic form in terms of either diameter or cross-sectional area as ( 1 ) + ( 2 ) = 1 (1) ( 1 ) 2 + ( 2 ) 2 = 1 (2) where the subscripts m, dd1 and dd2 denote the mother and daughter vessels respectively. τ is the scaling parameter which characterize the different scaling laws (13). The well-known Murray's law was the first scaling law proposed, where τ = 3 was derived by minimizing the combined viscous power dissipation and the metabolic power expenditure across the vascular network (40,41). Subsequent investigations (27,76) demonstrated that the main limitation of the Murray's law was that it considered each bifurcation in isolation instead of being part of a complete network and proposed successive improvements (based on the minimization of the cost of fluid conduction and fluid metabolism) culminating in the well-established HK law, with an exponent = 7 3 (23,24). In the analysis below, these two laws are considered in detail.

1D blood flow theory
The mathematical background of the 1D blood flow formulation has been extensively described in literature (14,36,37,58,59). Importantly, the forward and backward wave reflection coefficients at a bifurcation derive from the following system of conservation equations (mass and momentum) in three variables: cross-sectional area, pressure and velocity (A,p,v) where α is a non-dimensional correction factor for momentum flux, ρ is the blood density and κ represents the viscous resistance of the flow per unit length of vessel (60). This system is closed by a constitutive law relating pressure to area, derived from an elastic linear shell model (2), Here, A0 indicates the reference area and β the vessel material properties, which in turn are dependent on the Young's modulus E, the vessel thickness h and the Poisson's ratio ν which is usually taken as 0.5 since biological tissue is nearly incompressible (14,59).
The formulation outlined above for a single blood vessel can be then extended to a vascular network by imposing suitable coupling conditions at the vessel junctions (58,59). A common approach is to represent junctions as a single point and to disregard the impact of the branching angles and momentum loss, since it has been shown that they play only a minor role on wave propagation in the physiological range of pressure and velocity in the coronary vessels (14).

Wave transmission theory
In this work we employ a linearized wave transmission regime. Although this approximation is most accurate when the variation in the underlying hemodynamics is small, the tractability of the linear analysis can offer several important insights as demonstrated in Sherwin et al. (59). In this theory the reflection coefficient Rf at a bifurcation is defined as the ratio of the amplitude of the reflected and incident waves, where the characteristic impedance of a vessel Z0 relates the velocity or flow of a wave with the applied pressure and is defined as Consequently, the reflection coefficient Rf can be rewritten in terms of cross-sectional area and PWS as under the reasonable assumption of constant blood density. It is important to highlight that, from Eq. 7, the pulse wave speed in the reference configuration is Note that commonly in the definition of β (Eq. 6) the term A0 is included. However, since in this analysis A = A0 due to linearisation, it is preferable to introduce an area-independent measure β * . The backward reflection coefficients for the daughter vessels Rd1, Rd2 can be analogously derived by swapping in Eq. 10 the term related to the mother vessel 0,

Unified framework
To explore their intrinsic relationship, the scaling law formulation and the wave transmission theory have to be brought into a common framework. This can be achieved by recasting the mathematical formulation of the two frameworks in terms of the area ratio σ and the symmetry ratio γ (see Fig. 1), defined as Here onwards, σ and γ are used as the main parameters in our model derivation. Note that γ = 0 implies no bifurcation whereas γ = 1 describes a perfectly symmetrical bifurcation.

Figure 1 preferred placement
Starting from the generic form of a scaling law (Section 2.1) and considering that 1+ = 1 , Eq. 2 can be rewritten in terms of σ, γ as A closed form expression of σ in terms of γ can be then easily derived as Moving to the wave transmission theory (Section 2.3), assuming without loss of generality that A ≈ c ξ the forward reflection coefficient can be rewritten as after the variable substitution 1 − = 2 . It is then evident that imposing Rf = 0 in Eq. 17, thus implying no reflections for forward travelling waves, provides the theoretical σ − γ relationship of Eq. 16, which has been derived from the scaling law framework. This is a crucial point of the analysis and will be extensively discussed in the Result section. The backward reflection coefficients Rd1, Rd2 can be rewritten in terms of σ and γ following the same steps as presented above.

Impact of pulsatility
The derivation above has unified the framework of scaling laws with the wave transmission theory, and enables the link between branching structure (as characterized by a specific scaling law) and wave propagation  (28), the more general form of the characteristic impedance Z0, described in (3,44,46), should be considered when addressing the coronary circulation: where 0 ′ and 0 are functions of the Womersley number α (70,71). These functions have been introduced by Womersley (70) to obtain a more compact expression of the velocity profile in pulsatile flow (44). Full description and significance of each of these variables are briefly outlined in the Appendix, and more extensively in the references (44,70,71). Essentially, in this study the values of 0 ′ and 0 for different values of α are derived from those tabulated in Womersley (70).
Therefore, following the analysis presented in Brown (3), for a given Womersley number in the mother vessel αm, the corresponding values in the daughter vessels can be calculated as This reveals that at each junction α may vary heterogeneously from segment to segment as a function of the area ratio and symmetry ratio. The reflection coefficient Rf in Eq. 17 can be rewritten, following the same steps, to the casting resin for the human and canine heart. The injected casts were then allowed to harden for up to 24h at ambient temperature. Following this, the heart was immersed and the ventricular chambers were filled with carboxymethylcellulose sodium solvent (Brunschwig Chemie, The Netherlands) and Indian ink (Royal Talens, The Netherlands) mixture, and stored at −20 • C. The samples were then transferred to the cryomicrotome and sequentially sliced and imaged using fluorescent optical surface imaging (16,19). The image stack of each vasculature was then imported into our custom-developed automatic vascular extraction software and the geometrical representation graph was obtained. The three porcine hearts processed were resampled to isotropic voxel size of 58µm, 54µm and 53µm respectively, while the original voxel size were half these dimensions. The voxel size of the human heart and the canine heart were 58µm and 50µm respectively.

Figure 2 preferred placement
The main steps behind the custom 3D vessel segmentation pipeline are summarized in Figure 2. Note that the current work utilises an enhanced implementation, compared to the one presented in (16).
To improve the quality of the original image stack (4000x4000x2000 image slices) image restoration techniques were initially applied including deconvolution algorithms to tackle blurring (16). Subsequently, the range of vessel sizes typical of the coronary vasculature were detected separately by applying multiscale Hessian-based filters and then combined together (15). In addition, the multi-scale vesselness descriptor was expanded by complementing it with the addition of the Sato vesselness measure (53). The centerline extraction was then performed by applying the medial surface/axis thinning algorithm on the binarized volume (39).
Subsequently the vessel diameters were estimated and sub-voxel adjustment of the medial points were achieved by applying the Rayburst sampling algorithm (49) combined with a 3D core that follows the vessel centerline providing higher accuracy than the sphere fitting algorithm (16). The quality of the extracted vasculature network was further enhanced by a series of algorithms that prune out unphysiological structures across the vasculature, typically caused by cast leakage.

Validation of the unified framework
The newly developed theoretical framework presented above (see Section 2.4-2.5) is validated using anatomical measurements obtained from three high-resolution cryomicrotome volume images of ex-vivo porcine coronary vasculature, visualized in Fig. 3. These vascular trees provide anatomical measurements of ≈ 80000 bifurcations with a radius range of 1.6mm−50µm. The data were obtained by a custom-developed experimental setup and vascular extraction pipeline which were described in the previous section.

Figure 3 preferred placement
The high resolution cryomicrotome volume images provide sufficient anatomical detail to be analyzed both from a hierarchical point of view (radius of the mother vessel) as well as separating the territories perfused by the large epicardial vessels (LAD, LCx, RCA). In addition, for each of the porcine vascular networks (see Fig. 4), the myocardium (comprising both the left and right ventricle) has been manually segmented (using ImageJ (54)), and tetrahedralized (using CGAL (65)). Following this step, the Laplace equation has been solved with Dirichlet boundary conditions comprising 0 and 1 as boundary condition on the epicardial and endocardial surfaces to obtain a linear transmural gradient through the myocardial wall. The transmural depth was then assigned to each bifurcation using a kd-tree based search algorithm (to find the k nodes of the ventricular mesh 11 closest to a vascular node) with k = 10 thus offering an additional classification of the extracted coronary networks in terms of transmurality.

Figure 4 preferred placement
From this data, the statistical distributions of both σ and γ were initially calculated. The spatial distributions of the area and symmetry ratios ware then investigated to assess whether these quantities correlate with different vessel scales, perfused territories or transmural location. The σ-γ theoretical relationship (Eq. [16][17][18][19][20][21][22][23] for the different scaling laws considered (Murray's law and the HK law) was then compared with the measurements, including the impact of varying Womersley number. Following this analysis, the forward and backward reflection coefficients (Rf, Rd1 and Rd2) were calculated directly from the experimental data, using Eq.
10, and compared with the theoretical predictions. However, note that Eq. 10 depends on the knowledge of both the anatomical measurements and the PWS (or alternatively the material properties β * ) in each vessel.
Thus, the distribution of the PWS (or β * ) across the vasculature had to be hypothesized since it cannot be experimentally measured in every vessel at present. Therefore in the following, we examine several common hypotheses employed in the literature. Two simpler assumptions are to impose a constant uniform PWS (37) (which implies decreasing β * distally) or uniform material property (implying an increase in the PWS distally) where k1 = 2 × 10 7 g · s −2 · cm −1 , k2 = −22.53cm −1 and k3 = 8.65 × 10 5 g · s −2 · cm −1 . Each of these three approaches was applied on the extracted porcine coronary vasculatures to calculate Rf, Rd1 and Rd2 at each bifurcation.
Finally, the results derived from porcine vasculatures were compared with those from the other species (1 human and 1 canine vasculature) to investigate possible inter-species variation.

Unified framework
The theoretical analysis described in the previous section (Eq. 14-17) has multiple implications. It demonstrates for the first time, to our knowledge, the equivalence between the adherence of a vascular network to a particular scaling law and the well-matchedness of its bifurcations, thereby closing the loop between vascular design and wave propagation dynamics. Moreover, if well-matchedness is assumed, then a σ-γ closed form analytical relationship is derived (Eq. 16) and can be used for validation, since it depends solely on the anatomical measurements without requiring the knowledge of the pulse wave speed or the material properties in each segment, which are not experimentally measurable currently.

Figure 5 preferred placement
In Fig. 5a the theoretical σ-γ relationships for the two different scaling laws considered are visualized, along with experimental data points measured from human epicardial coronary vessels. If well-matchedness is assumed a priori then both Murray's law and HK law predict that more symmetrical bifurcations (γ > 0.5) would exhibit higher area ratios -this can be tested against experimental measurements. It can be seen that Murray's law predicts higher area ratios for a given symmetry ratio compared to the HK law, and neither law captures the spread of the data points adequately. Furthermore, the measurements from Finet et al. (12) and Russell et al. (52) cover a relatively small area ratio range (γ ≈ 0.8) and are limited to large vessels (d > 2 mm), which motivates the need for further experimental investigation.
The impact of Womersley number α on the σ-γ relationship is visualized for the HK law (Fig. 5c) and Murray's law (Fig. 5d)

Vascular extraction results
Before moving into the data analysis, the reliability of the vasculature extraction pipeline is compared against previous literature here. The percentage of bifurcations and n-furcations (with n > 2) in each extracted porcine vasculatures is ≈ 95% and ≈ 4% respectively, which is in good agreement with Kassab et al. (29). The trifurcations are excluded in the subsequent analysis due to their limited numbers in the coronary vasculature (29,67). Moreover, as visualized in Fig. 6, the relationship between the mother and daughter vessels diameters in each network compare well with that of VanBavel and Spaan (67) demonstrating the reliability of the radius estimation step in the vascular extraction pipeline. Finally, since the statistical distributions of the area ratio and symmetry ratio visualized in Fig. 6 (last two columns), agree well with VanBavel and Spaan (67) and are qualitatively consistent between samples, the datasets were merged together for the remainder of the analysis. There is a statistically significant (p<0.01, Wilcoxon signed-rank test) increase in the area ratio (from 1.05 to 1.4) and symmetry ratio (from 0.3 to 0.7) moving from the large epicardial vessels to the small arteries and arterioles, as shown in Fig. 7. Furthermore, the percentage of asymmetric bifurcation (γ < 0.5) strongly decreases from ≈ 80% in the large vessels to ≈ 15% in the smaller ones ( Fig. 7c-d). However, note that there is no statistically significant difference in area ratio and symmetry ratio when comparing different perfused territories (Fig. 7a-c) or transmural location (Fig. 7d-f) within the same vessel size range. These findings remain unchanged when the Rm interval divisions are adjusted, as changing the mother vessel radius ranges up to 20% causes only a minimal variation (< 5%) in the metrics analysed.

Validation of theory
The scaling law coefficient τ was fitted, using the spatial distribution of σ and γ across the coronary vasculature (Eq. 16). The fitting process was repeated 50 times for each cluster region and the solution τ providing the minimum mean squared error was chosen. As clearly visible in Fig. 8a, there is a statistically significant (p < 0.01) increase in the scaling law coefficient τ moving from large (τ ≈ 2.25) to small (τ ≈ 3.5)  To summarize, when the pulsatility of flow was disregarded in the analysis of coronary branching patterns, the analysis concluded that different scaling laws apply to different vessel scales (see Fig. 7). However, when the theoretical framework of scaling law was extended to incorporate the Womersley number, the HK law satisfactorily agreed with the anatomical data for vessels of radius larger than 0.1mm. The same conclusions were obtained when the vasculature was clustered according to the transmural depths, although the details are not shown here.

Reflection coefficients across the vasculature
The forward (Rf) and backward (Rd1, Rd2) reflection coefficients were calculated for each mother vessel radius, perfusion territories or transmural location, for the three parametrizing approaches described in Section 2.7 (uniform PWS, uniform β * or a experimentally-fitted relationship). The outcomes were somewhat surprising in that the proximal Rf, Rd1 and Rd2 were largely similar (< 10% variation) for all parametrization scenarios despite the fact that the assumed distribution in the PWS/β * varied significantly. Moreover, similar results were observed when clustering the bifurcations in terms of the mother vessel radius (Rm) or the Weibel generation number (29). As shown in Fig. 10a-d for the constant β * assumption, the forward reflection coefficient is small (Rf ≈ 0) thus supporting once more the well-matchedness hypothesis of the coronary bifurcations for vessels with a radius Rm > 0.1. Smaller vessels have a slightly more negative forward reflection coefficient of approximately -0.1. Moreover, due to the observed increase in symmetry ratio and area ratio for smaller vessels, the backward reflection coefficient Rd1 of the larger daughter vessel becomes more negative while the smaller daughter vessel coefficient becomes less negative. This is in good agreement with the theoretical predictions visualized in Fig. 5.

Species Comparison
The analysis performed above was then repeated for the canine and human vasculatures (Fig. 11), to assess possible inter-species variation. Although the superficial appearance suggested a potentially different relationship (the human coronary vasculature is more tortuous compared to the canine and porcine counterparts), the same qualitative trends and quantitative results were confirmed, both regarding the σ, γ distribution across the vasculature and the satisfactory prediction of the σ − γ relationship by the HK law, when the impact of pulsatility was taken into account.

Discussion
The key outcomes of the current study are as follows -the augmentation of the scaling laws with the wave transmission theory enabled the implications of specific scaling laws to be studied in an integrated manner, in terms of bifurcation morphology and wave reflections. The consequent validation using the pulsatile regime model offered a plausible explanation for the heterogeneity of branching morphology that was previously reported in literature (12,52,77). More importantly, the theoretical results and experimental data supported the hypothesis that the forward well-matchedness is a salient feature of coronary vascular structure.
Furthermore, this observation has been shown to be robust to different assumed distributions of vessel material properties and PWS that have been trialed. This in turn implied that the backward well-matchedness is selective; that is, transmission is only preferred from the larger of the daughter vessels, and suppressed for the waves travelling up from the smaller daughter vessel. In the discussions below, we examine the theoretical and experimental aspects of the current work and the physiological implications of our findings in greater detail.

What new findings does the proposed framework offer?
In this work, the unification of scaling laws and wave transmission theories led to a demonstration that assessing the well-matchedness of the vessel branching pattern is equivalent to assessing the validity of a specific scaling law. It was shown that recasting the mathematical formulations of both frameworks in terms of area ratio σ and symmetry ratio γ (see Section 2.4) provides the same σ-γ relationship (Eq. 16) under the assumption of zero reflection coefficient for forward travelling waves. This is crucial for two main reasons. For validation purposes, two main predictions can be deduced from the newly introduced theoretical framework and compared with the experimental measurements. Firstly, more symmetrical bifurcations (larger γ) have higher area ratios (larger σ), as visualized in Fig. 5a. Secondly, smaller vessels have higher area ratios for a given fixed symmetry ratio consistent with the expected reduction in Womersley number (see Fig.   5c-d). Both predictions are confirmed by the high-resolution measurements (of around 80000 bifurcations with a 40-60µm minimum vessel diameter, see Fig. 9).

Are the existing scaling laws valid?
The analysis excluding the influence of pulsatility has called into question the validity of a single scaling law applied to the whole vessel network. When fitting the scaling coefficient τ, different values for different vessel sizes were obtained (see Fig. 8), suggesting that scaling relationship itself may depend on the vessel diameters.
This dependence of the scaling law parameter τ on the diameter of the vessels has been experimentally investigated and reported in literature. Specifically, in VanBavel and Spaan (67), τ was found to be 2.82, 2.5 and 2.35 for a mother diameter of < 40µm, 40µm − 200µm and > 200µm, respectively. Moreover, the decrease of τ for an increase in the vessel radius (moving from small to large vessel) has been confirmed by other studies (1,67). These results compare well with our findings (Fig. 8).
On the other hand, the full analysis including pulsatility demonstrated the possibility that the HK law is compatible with the measurements when Rm > 0.1, if one takes into account the range of Womersley numbers found in vivo. This further implies well-matchedness throughout the network, according to our analysis. This is because the Womersley regime acknowledges the heterogeneous variation of physiological σ − γ relationships (Fig. 5) dependent on the local hemodynamic conditions, rather than demanding a single curve corresponding to non-oscillatory conditions. However, note that we have only examined the branching pattern aspect of the scaling law. Similar variability in the scaling relationships regarding other morphometric parameters (tissue volume-vascular volume relationship) have been recently examined in van Horssen et al.
(20), and remains a task for future exploration.
To summarise, even when pulsatility was taken into account, the HK law satisfactorily fitted the data and this implies well-matchedness for coronary bifurcations. However, it is important to underline that the conclusions drawn so far are valid down to a vessel scale of Rm ≈ 0.1 mm. In fact, as can be seen in Fig. 9, for smaller vessels (particularly the asymmetrical ones) the measurements do not fit the theoretical predictions.
This may be due to the fact that the anatomical reconstruction at this scale is subject to greater error since the vessel dimension approaches the image resolution limit. However it is also likely that wave propagation, on which the newly developed theoretical framework is based, plays only a minor role at this scale. In fact, the vessel segments with Rm < 0.1 mm have an average length of approximately 14mm, and a PWS derived from our analysis of approximately 20 − 25 m/s (depending on the parametrization approach assumed). This implies that a wave would travel through a vessel in less than 1ms (not measurable by existing techniques and which would require a sampling frequency >1 KHz). In this regime, the previously proposed Windkessel approximation (i.e. infinite PWS) of the coronary vasculature may be more appropriate (30).
Finally, it is important to recall once more that each proposed scaling law implies a specific branching pattern, that has never been evaluated in terms of its consequence to wave propagation. Specifically, once a certain scaling law is assumed, the consequent forward and backward reflection coefficients can be calculated at each bifurcation thus enabling the wave reflection pattern to be assessed.

What does the current evidence show regarding coronary wave reflection properties?
The validation of the newly introduced theoretical framework supports the hypothesis that the coronary bifurcations are well-matched for forward travelling waves, whereas the backward travelling waves are significantly damped, increasingly so at smaller scales (R<0.3mm). This conclusion can be derived firstly by the observation, in Fig. 9, that the anatomical measurements in vessels with Rm > 0.1mm agree with the theoretical predictions that have been obtained under the assumption of well-matchedness. Secondly, the reflection coefficient for forward travelling waves has been found to be Rf ≈ 0 for the analysed scales, regardless of the parametrization approach used (constant material properties, constant PWS or experimental formula), with no significant differences in terms of transmurality or perfused territories (see Fig. 10).

18
The well-matchedness for forward travelling waves subsequently implies ill-matchedness (strongly negative reflection coefficients Rd1,d2) for at least one of the backward travelling waves at the junction, especially for more asymmetrical bifurcations (Fig. 10).
It is noteworthy that the distribution of the forward and backward reflection coefficients has been calculated using the linearized formula derived from Eq. 17, instead of Eq. 23, which takes into account pulsatility. This choice was motivated by the need to minimize the number of assumptions required for calculating the reflection coefficient at each bifurcation (see Section 2.7). In fact, Eq. 17 assumes a priori the distribution of PWS (or material properties β * ) across the vasculature. Using Eq. 23 would additionally require the assumption of Womersley number distribution across the network (or to derive it from previous publications (28)). The fact that different PWS distributions led to the identical outcome of well-matchedness for forward travelling waves strongly suggests that this is a salient feature of the coronary vasculature.
Furthermore, this conclusion is in good agreement with Huo and Kassab (22) model of pulsatile blood flow in the entire coronary tree, in which it was observed that the damping of the pressure waveform is rather small for larger vessels and becomes more pronounced for vessels smaller than 100µm. However, neither the reflection coefficient distribution nor the PWS distribution across the vasculature have been reported, thus a direct comparison with our current results is not possible.

How are these results relevant to clinical challenges?
Waveforms in the coronary vessel have been studied for their clinical diagnostic potential (62,63). Many recent investigations have employed techniques such as wave intensity analysis (47), that permits specific features of the wave to be identified and ascribed to the spatio-temporal sequence of the coupled cardiaccoronary cycle and their deviation from the norm under diseased conditions. In particular, the role of the backward expansion wave (BEW) has received attention due to its alteration in hypertrophic cardiomyopathy (8), aortic stenosis (7), heart failure (34) and prognosis following myocardial infarction (9) from which novel insights regarding the specific pathology have been derived.
However, while the origin of the forward travelling waves is reliably ascribable to the ejection/suction of the LV, the mechanisms underlying the backward travelling waves are yet to be fully described thus limiting the conclusions derived by cWIA. Specifically, in the early clinical research studies, the origin of the backward travelling waves has been putatively ascribed at the microcirculatory level, following the relief of systolic myocardial compressions (7,8,17,34,55). However, there is no specific evidence quantifying at which depth of the network such relaxation is taking place to give rise to the detected waveforms -in other words, the location of the horizon of observable backward waves remains unidentified.
A more recent study, investigating the variation of the cWIA profile during the Valsalva maneuver led to the conclusion that coronary wave energy is directly influenced by cardiac mechanical factors (maximal and minimal LV pressure rate) and poorly related to mean coronary flow (51), thus suggesting the compression/decompression of the intramural vessel as the main mechanism driving the backward travelling waves. Our current estimates of the backward reflection coefficients (Fig. 10) not only support this finding (since waves generated at the microcirculatory level would hardly reach the inlet of the epicardial vessels due to the strongly negative reflection coefficients Rd1,d2) but additionally indicate that the observable wave horizon may be markedly proximal in the coronary network. For example, at Rd1 = −0.3, only 70% of the waves from the larger daughter vessel would be transmitted up to the parent at each junction. The smaller daughter would effectively contribute no detectable waves past one or two generations.
These observations motivate a reexamination of the existing physiological and pathophysiological conclusions derived by applying cWIA. For instance, the lack of diastolic dominance of flow velocity in the right coronary artery compared to the left main stem, associated with a less prominent BEW, is more likely due to the smaller compressive force on the intramural vessels of the right ventricle than on the coronary microvasculature, as originally suggested (17). More importantly, the prognostic value of the BEW in predicting myocardial recovery post-infarction may derive from its direct relationship with the LV contractility more than its capability to assess the integrity of the microcirculation, post-MI (9).
In summary, the backward travelling waves are thought to be generated by the In addition, the well-matchedness of the coronary bifurcations in the large epicardial vessels is an important feature that must be taken into consideration when designing and implanting a stent or a coronary bypass graft (13). The diameter and the material properties should be carefully chosen to maintain the wave transmission properties, otherwise a significant rise in wave reflection at the top of the coronary vasculature may lead to increased shear stress burden and resistance to flow.

Relevance for one-dimensional modelling of coronary blood flow
Though one-dimensional blood flow models have gained increasing importance in modelling applications due to their efficiency and the ability to accurately represent wave propagation phenomena (38,66), the parametrization of the vascular network remains an unsolved challenge. Overcome by the difficulties of directly measuring the parameters, many studies adopt the strategy of manually tuning individual segments to achieve idealized wave behavior. Our proposed theory can be used to assess different approaches for assigning wall stiffness parameters to each segment of the network in simulation studies.
For instance, in the low α regime, imposing a uniform PWS or material parameter β * implies τ = 2 and τ = 2.5 respectively (3). Substituting these values in Eq. 16 provides the corresponding σ − γ curves. In qualitative terms this would result, for the uniform PWS case, in a constant value of the area ratio σ for all the possible symmetry ratio values, which is not supported by the experimental measurements. On the contrary, the uniform β * assumption would correctly predict increasing area ratios for higher symmetry ratios.
Although introducing the pulsatility (Eq. 23) to this analysis bestows correct behavior to the uniform PWS scheme, the uniform β * assumption still leads to a superior fit to the experimental data. Despite these differences, all three approaches predict a similar PWS and β * distribution for vessels of radius larger than approximately 0.5mm. For vessels smaller than this scale, the experimental formula and the uniform β * hypothesis predict a strong increase in the PWS in contrast to the uniform PWS hypothesis which predicts a decrease in β * across the generations.
What are the major sources of experimental error?
Significant efforts have been made to ensure the reliability of the experimental procedure and vascular extraction pipeline. In our experiments, a strict protocol was followed to minimize the variability in the data acquisition, and repeatability was demonstrated by the consistent outcomes between the three porcine samples. However, it is possible that adenosine may have had a non-negligible systematic impact on the anatomical measurements. Adenosine, as a vasodilator, may cause a dose-dependent dilation in the coronary vessels beyond a physiological range, leading to a significant drop in resistance (11,32,33). Nevertheless, its effects are strongest in the small intramyocardial vessels (R< 50µm) (11), below the scale where wave propagation is likely to play a strong role as previously discussed. Most importantly, our analysis is based on the area and symmetry ratios, which depend on diameter measurements in close proximity (adjacent to each bifurcation). Thus, it is unlikely that adenosine would have caused a significant bias in our measurements.
As described in van Horssen et al. (20), vascular reconstruction from cryomicrotome imaging is limited by the degree of penetration of the casting material and by the imaging resolution of the experimental setup.
However, previous validation studies have shown that the casting material reliably penetrates vessels down to 10µm, which is far below the resolution required for this study.

Limitations of the study
In a previous investigation it was highlighted that even large changes in the scaling law parameter (τ) led to only a small impact on the cost function of the energy minimization (57). This curtailed the possibility of distinguishing between different scaling laws using experimental data. On the contrary, in our model changing τ within the range of interest (from 7 3 to 3) caused, for a given symmetry ratio, a substantial variation in the theoretically predicted area ratio (see Fig. 5 and Fig. 9) thus avoiding the ambiguities of the previous analyses.
The validity of applying the Womersley analysis in the coronary vasculature, on which this work is based, is a matter of possible controversy. In fact, it can be argued that the distance between subsequent bifurcations in the coronary vasculature is not sufficient to achieve the implied (Womersley) velocity profile. The question of the optimal theoretical approach to describe the velocity profile in the coronary flow is still under debate in the scientific community as remarked in van de Vosse and Stergiopulos (66), and is hampered by the complexity of measuring it in in vivo settings. Nevertheless, multiple three-dimensional in vitro and in silico studies in the literature illustrated the benefit of employing the idealized Womersley profile to better describe local hemodynamic conditions and wall shear stress distributions in the large epicardial vessels (18,31). This was corroborated by our study which demonstrated a better agreement between the theoretical predictions and the experimental measurements (see Fig. 9) by adopting the Womersley analysis.
The analysis in the present work was mainly based around diastolic conditions, consistent with the experimental set up from which the anatomical data were collected. The rationale is that for consideration of general structural design of a vascular network, the first order effects are to be found in the dominant operating condition i.e. diastole, during which around 80% of the arterial flow occurs. The same logic can be applied to variation of heart rate from resting norm (which may alter the pulsatility) and systolic-diastolic variations (resulting in anatomical changes as well as pressure difference, leading to altered pulse wave speed in segments). Nevertheless, it would be instructive as a future exploration, to examine how robust the current findings will remain when the cardiovascular system is subjected to varying loads.
Finally, in our work the effects of vessel tapering and non-linearities on wave propagation have been ignored. It is known that non-linearity in the pressure-area relationship affects wave-propagation and that tapering in the geometrical or material properties generates wave reflections. However, our preliminary numerical simulations (not shown) in single tapering vessels in approximate dimensions of a large epicardial vessel demonstrated that tapering and non-linearities counterbalance each other and only marginally impact the forward and backward travelling wave magnitudes (< 5%).

Conclusion
A new theoretical framework which unifies the wave transmission theory with the scaling law formulation has been introduced and validated for the first time. Including pulsatility in the framework broadened the range of measured bifurcation anatomy that can be explained, compared to the existing models. 22 The validation demonstrated that the branching structure in the healthy coronary vasculature is wellmatched thus favoring forward wave propagation. This feature holds robustly in the proximal network in several species (porcine, canine and human), regardless of the various assumptions posed on the distribution of pulse wave speeds (equivalently, the wall stiffness) among the segments. This showed, conversely, that backward travelling waves are strongly damped when passing through a bifurcation thus indicating that the horizon of observable backward waves could be significantly more proximal than previously thought in the coronary vasculature.
In addition, new data providing branching patterns across the coronary vasculature has been described in a high level of detail showing that smaller vessels have higher area ratios, higher symmetry ratios and a higher numbers of asymmetrical bifurcations than the large epicardial ones. However, no significant differences was observed in terms of transmurality or perfused territory (LAD, LCx or RCA).

Appendix
The functions 0 ′ and 0 have been introduced by Womersley (70) to obtain a more compact expression of the velocity profile in pulsatile flow (44). More specifically, this involves the substitutions        The σ-γ relationship for the two scaling laws considered is visualized. Murray's law predicts higher area ratios for a given symmetry ratio than the HK law. (b) The reflection coefficients for backward travelling waves for both the larger and smaller daughter vessel are visualized. It is evident how asymmetrical bifurcations (γ < 0.5) strongly disadvantage backward wave propagation (highly negative reflection coefficient) in the smaller daughter vessel. (c-d) The shaded region shows the impact of varying Womersley number on the σ-γ relationship (Eq. 23) for both the HK law (c) and Murray's law (d). α > 10 represents the case when the effect of pulsatility is disregarded (Eq. 16). For lower α the theory predicts that the well-matchedness hypothesis is fulfilled at higher area ratios for a given symmetry ratio.

Figure 6:
For each porcine vasculatures analysed the relationships between the mother diameter (dm) and the large (dd1) and small (dd2) daughter vessel diameter are visualized using a log-log scale in the first and second column respectively. The trend is in good agreement with the literature (65). The histograms of the area ratio σ and symmetry ratio γ are presented in the third and fourth column. The distribution is qualitatively consistent between vasculatures.

Figure 7:
The mean and standard error of the area ratio σ (a-d), symmetry ratio γ (b-e) and percentage of asymmetrical bifurcations (c-f) are visualized for each vessel cluster. The vasculature has been classified in terms of mother vessel radius, perfused territory or transmural location. It is clear that, as the vessel radius decreases all these three variables increase. Interestingly, no difference was found when comparing the same vessel size in different perfused territories or transmural locations.

Figure 8:
As the mother radius diameter decreases, the scaling law coefficient τ (see Eq. 2) increases. This estimation of τ acquired under the assumption of steady flow regime, implies that different scaling law coefficients may pertain to different vascular regions. Figure 9: The area ratio-symmetry ratio relationship for the Womersley number range typical of the coronary vasculature is visualized for the HK law (shaded grey) and Murray's law (striped black) overlapped, for the different mother vessel radii (Rm) intervals considered. Firstly, the theoretical prediction of higher area ratios for higher symmetry ratios is confirmed by the measurements, up to a scale of Rm > 0.1 where the data exhibit a plateau of σ ≈ 1.4. When the varying degree of pulsatility is taken into account, the HK law satisfactorily agree with the data for the large, medium and small vessel (a-c) range considered. Moreover, it predicts an increase in the σ − γ relationship due to a decreasing Womersley number, which is a consequence of decreasing vessel radius. The theoretical predictions for Murray's law (striped black area) do not match the measurements for the large and medium scales.

Figure 10:
The wave reflection coefficient for the forward Rf and backward travelling waves Rd1,d2 were calculated for the uniform β * hypothesis clustering the vasculature in terms of mother vessel radius Rm (a-c) and Weibel generation number (d-f). Note that the Weibel generation number starts from 1 at the first bifurcation and increases by 1 at each consequent bifurcation. Similar results were obtained using the hypotheses of uniform PWS as well as the experimentally fitted formula (45). Note that in the rage of Rm > 0.1 the forward reflection coefficient is approximately zero and supports the well-matchedness hypothesis of the coronary bifurcations. Moreover, smaller vessels were found to exhibit less negative reflection coefficient in the smaller daughter vessel Rd2 at the expenses of the larger daughter for which Rd1 became more extreme.

Figure 11:
The extracted canine and human coronary vasculature are visualized. The LAD, LCx, and RCA perfused territory are highlighted in red, gold, and silver respectively. Note that the human vasculature is more tortuous than the porcine one.