Cohort-based multiscale analysis of hemodynamic-driven growth and remodeling of the embryonic pharyngeal arch arteries

Congenital heart defects (CHDs) constitute the most severe form of a birth defect, accounting for 24% of deaths from developmental abnormalities (Mozaffarian et al., 2016). Abnormal pharyngeal arch artery (PAA) morphogenesis is associated with over 50% of clinically presented congenital heart defects (Go et al., 2013). Hemodynamics and proper flow patterns during cardiac development have been acknowledged as major drivers of vascular growth and remodeling (Lin and Taber, 1995; Jaffee, 1965; Taber and Eggers, 1996). Blood flow and its associated forces determine the mechanical environment of a subject, which in turn has implications at the cellular and molecular level. Wall shear stress (WSS) is sensed by endothelial cells and has been shown to regulate vascular remodeling (Brownlee and Langille, 1991; Kamiya and Togawa, 1980; Langille and Donnell, 1986), though exact thresholds for such changes are likely regional and stage specific. Pressure has been shown to modulate vessel shape and stretch (Banerjee et al., 2015; Taber, 1995; Li et al., 2004). Both mechanical and epigenetic perturbations in flow have been implicated in various congenital heart defects (Hu et al., 2009; Sedmera et al., 1999; Yashiro et al., 2007), including aberrant development of the great vessels (Abu-issa et al., 2002; Hogers et al., 1997, 1999; Rychter, 1962; Wendling et al., 2000). Yet despite all these associative relationships, the numeric relationship between mechanical forces and abnormal morphology, particularly at embryonic stages, has remained elusive. The difficulty of this quantification can be attributed to a lack of cohort-based analysis in numerical studies designed to quantify in vivo findings in which ranges for norms are reported. In order to render early CHD risk assessment and intervention possible, there is a need to establish cohort-based variations in flow and morphology. Only after the establishment of such a range will we as a field be able to assess its connection to malformation risk.

To date, quantitative studies of PAA development rather than qualitative surveys of PAA morphology are few (Hu et al., 2009; Kowalski et al., 2014; Leatherbury et al., 1990a,b; Lindsey et al., 2015; Wang et al., 2009). Leatherbury et al. (1990a,b) assessed PAA function through microcinephotography images and pressure and diameter measurements at the level of the ventricle and dorsal aorta (DoA), coupled with qualitative observations of the arches themselves. Hu et al. (2009) reported HH24 pulsatile Doppler flow measurements for PAA III and IV pairs, but were unable to measure the PAA VI pair owing to its position in the pharyngeal mesenchyme. No other published reports of experimentally measured PAA individual arch artery velocities exist. Similarly, in vivo pressure measurements have been limited to the DoA (Yoshigi et al., 2000; Zahka et al., 1989). Computational examinations of the arch arteries (Kowalski et al., 2012, 2013; Lindsey et al., 2015; Wang et al., 2009) have quantified blood flow and spatial variations in WSS across the entire PAA system but have not reported pressure maps. Based on idealized composite models, those studies elucidated the global role of flow and WSS in embryonic heart development but not local subject-specific hemodynamics that could reveal potential zones of high or low variability. The PAA network with its three lateral arch pairs serves as a highly informative model for quantitative variability studies, both locally, via segment by segment quantification across arches, and globally, via comparisons between embryos.

Before formation of the pulmonary artery and aortic arch, blood exits the heart through six bilaterally paired pharyngeal arch arteries. Growth and remodeling of these vessels consist of a complex sequence of events that involves the emergence, disappearance, elongation, rotation and twisting of the arch pairs, as they remodel into their mature adult forms. In the chick embryo, the process of growth to maturation of the arch pairs largely occurs over the course of 7 days. The mature pulmonary artery and aortic arch are present by HH36 (day 10). Throughout maturation of the arch arteries, three of the six arch pairs persist: PAA III, PAA IV and PAA VI. These final arch pairs are present by HH24 (day 4), at which point the outflow tract has shifted to its definitive position, ventral to the right atrium (Manner, 2000). PAA III contributes to the brachiocephalic and common carotid arteries. The right lateral PAA IV contributes to the adult aortic arch. The left lateral PAA IV regresses. PAA VI contributes to the ductus arteriosus and parts of the central pulmonary arteries. Of the transitory arch pairs, PAA I and II remodel into capillary beds, whereas PAA V is fleetingly seen as a segment of PAA VI (Hiruma and Hirakow, 1995; Hiruma et al., 2002). This sequence of growth and remodeling is conserved across species, although in mammals, the dominant side is reversed, with the left lateral PAA IV contributing to the mature aortic arch and the right lateral PAA IV forming part of the subclavian artery. Perturbations to this programmed sequence of events, can lead to severe cardiac abnormalities.

In this present study, we have developed a novel imaging and multiscale computational model pipeline to quantify and analyze local, segment by segment, stage-specific and subject-specific PAA hemodynamic and morphological changes during a critical window of development at HH18, HH24 and HH26 (Fig. 1), looking at the natural variation that may present itself at a particular stage. WSS and pressure are quantified in the different subjects and across stages, along with area, length and other morphological parameters, to determine which drive growth. We perform for the first time in the PAAs, a group-based parallel numerical study, using detailed segmentations of stage-specific arch artery systems obtained from in vivo measurements. We offer a new way of analyzing morphogenesis through the lens of quantitative variability analysis. Ultimately, establishing norms for critical windows of development that take into account stage-specific variability will help identify the origins of congenital malformations and provide insights into the optimal timing for restorative interventions.

In order to assess hemodynamic-driven growth, both morphological and hemodynamic parameters must be thoroughly categorized and a relationship between the two established. Characterizing morphological and hemodynamic growth patterns requires a method to track changes in PAA morphology and blood flow dynamics over a period of time. The more local the tracking, the less assumptions need to be made and the more variability can be analyzed. Our ‘high-resolution imaging and flow measurement to multiscale numerical analysis pipeline’ (Fig. 2) consists of two separate subject-specific numerical models obtained from 3D reconstructions of in vivo morphology. The 3D analysis allows for quantitative measurements of qualitative or previously unchecked observations. To create the pipeline inputs, embryos were open cultured and imaged via ultrasound. Pulsed-wave Doppler velocity curves were gathered for HH18, HH24 and HH26 embryos through B-mode guided videos. Small cohorts (n=5 per cohort) of embryos from each stage were perfusion fixed, stained with a potassium iodide and iodine solution and imaged by way of 3-4 µm voxel nano-CT scans. Resulting high-resolution image stacks were reconstructed and used as a basis for the numerical models. Doppler curves were then used as inputs into computational fluid dynamic (CFD) models to map out flow dynamics along the PAA manifold. The mapping of flow dynamics is important as blood flow and its resulting forces induce regional changes on the cellular and molecular level that are important for vascular morphogenesis. To make use of the detailed geometrical and hemodynamic parameters along the arch arteries of the 3D models, each arch was divided into ten equally spaced individual sections. This was achieved by establishing centerlines and perpendicular cross-sections along arch arteries and calculating desired parameters at each point. The division of arch arteries into segments allowed for the tracking of local growth within groups and across stages by enabling the comparison of growth between relative positions across multiple sources.

PAA evolution was first evaluated from HH18 to HH24, as it marks a critical remodeling period, notably the growing in of PAA IV, the primitive segment of the aortic arch, and the regression of PAA II into capillary beds. From HH18 to HH24, the outflow tract (OFT) itself shifts in position, moving ventral to the right atrium, and loses its once tubular character by HH24 (Manner, 2000). The embryo is rapidly growing and must account for a 108% change in blood volume (Rychter et al., 1955). In order for the PAAs to accommodate a larger blood volume, the arch arteries themselves must grow in length, diameter or cross-sectional area (CSA). Detailing these changes, allows us to understand how this growth comes about. Through the 3D reconstructions of in vivo morphology, calculation of curved arch length taken along the arch artery centerline was possible (Fig. 3). At HH18, PAA II, which will dissolve into capillary beds, is by far the longest arch, at 0.17±0.04 cm (right) and 0.17±0.03 cm (left), longer than its HH24 PAA III and IV (right and left) cranial counterparts and just under the PAA VIR and VIL (VI right and left, respectively) counterparts at HH24. The 3D reconstructions also revealed stage-specific shape and CSA changes along the arch arteries (Fig. 3A,B). Among HH18 embryos, vessel CSA are elongated along the horizontal (R↔L, Fig. 1) axis. By HH24, arch artery cross-sections become more round and less elongated in the horizontal direction. Centerline calculations allowed for evaluation of the way vessels branch off the aortic sac and its evolution over time. The angle at which a vessel branches off its parent vessel influences the immediate flow properties of the fluid moving through the conduits (Sinibaldi and Romano, 2017). This angle also dictates the amount of force (in the form of pressure) necessary for flow to change direction and embark on a new fluid path. In general, a smaller angle favors fewer secondary flows and requires a smaller pressure gradient to drive the flow through the junction. A branch-splitting angle (Fung, 1997) was calculated for each of the six arch arteries based upon the centerlines of the vessel of origin and the centerline of the new branch. Fig. 3D displays the branch splitting angles for the split of the DoA into a left and right branch, as well as well as the OFT inlet split into the six PAA vessels. For PAA III and VI (right and left), the branch splitting angle decreases with growth, whereas the angle of PAA IV (right and left) increases from HH18 to HH24. The angle can be thought of as the spreading or initial spacing between the vessels in the cranial caudal direction. As the vessels mature, they become more parallel, more in-line with each other and closer together, rather than splayed in the cranial-caudal direction. DoA branch splitting increases between HH18 and HH24. The left branch angle remains slightly more elevated than the right branch angle.

Fluid dynamic parameters are ultimately dependent on flow velocity and flow distribution throughout the arteries. To understand how flow dynamics shift throughout development, we examined changes in flow distribution over time/stages (Fig. 4) through the 3D numerical simulation. For HH18, the right and left third arch consistently receive the largest amount of flow in the developing HH18 embryo, receiving between 79% and 92% of flow through the arches over the course of one cardiac cycle. The right side receives more flow in two of the five HH18 embryos, whereas the left side received an average of 15% more flow in the three other cases (Table S1A). Within the HH24 cohort, the embryo obtains a more balanced flow distribution across the arches. Rather than one arch pair uniformly receiving upwards of 73% of the flow, HH24 embryos have an average of 42% (±12.9%), 39.3% (±4.1%) and 18.5% (±13.2%) mean flow distribution across the arch pairs (III, IV and VI, respectively). In HH24 embryos, the right fourth arch receives over 25% of the flow in three out of five embryos (Table S1B). Surprisingly, the left fourth arch received over 25% of the flow in the remaining two embryos and flow dominance was shifted to the left side in these embryos. When comparing arch pairs across stages, relative standard deviation (RSD) (standard deviation normalized to the mean value) stands at 62.9%, 6.2% and 35.7% for HH18 PAA II, III and IV pairs, respectively, and at 30.5%, 10.4% and 73% for HH24 PAA III, IV and VI pairs. These data suggest that more-malleable HH18 II and IV arches, in terms of flow distribution, permit a more controlled PAA III pair, whereas more malleable HH24 III and VI arch pairs permit a more controlled PAA IV pair.

Hemodynamic-driven (or mediated) growth could arise from two different mechanisms. Transmural pressure or WSS could lead to increases or decreases in vessel area or diameter. Using the 3D numerical simulation, we were able to analyze pressure and WSS values across PAAs (Fig. S1). HH18 peak pressure values vary across the cohort of five (Fig. S1A), but pressure patterns remain consistent: a peak that originates in the OFT inlet and dissipates cranially up the aortic sac and laterally through the arches, with the most caudal arch pair (IVR and IVL) absorbing the highest pressure per arch length (Fig. 4B). Peak pressure varies between 5100 dyne cm⁻² (3.83 mmHg) and 2900 dyne cm⁻² (2.17 mmHg) at peak flow (Fig. 4B), with an average arch pressure of 3387 (±128) dyne cm⁻². HH24 embryos mark a rapid increase in pressure, but conserve the same pressure dissipation pattern, with the most caudal arch pair absorbing the largest amount of pressure per arch length (PAA VIR and VIL) (Fig. S1B). HH24 pressure magnitude varies between 7500 dyne cm⁻² (5.63 mmHg) and 4650 dyne cm⁻² (3.49 mmHg) at peak flow (Fig. 4B). Average arch pressure increases to 6072 (±389) dyne cm⁻², an increase of 79% from the HH18 cohort. HH18 peak WSS is 48.15 (±23.59) dyne cm⁻² across embryos. Spikes in WSS can be seen at the inlet in the aortic sac near the third PAA entrance, as well as laterally along the length of the arch (Fig. 4B, Fig. S1D). HH24 embryos exhibit an increase in spatial frequency of WSS spikes with a per arch WSS magnitude of 94.74 (±5.58 dyne cm⁻²) at peak flow (Fig. 4B, Fig. S1E). Although the per arch magnitude of WSS greatly increased, HH24 WSS magnitude remains on the same scale as that of HH18 embryos.

To distinguish between pressure and WSS as drivers of PAA remodeling, we quantified morphology changes from CFD-aided 3D models using data obtained from the point-to-point segment markers (n=5 embryos per stage, 10 markers per arch, 6 arches per embryo). Statistically significant changes between the relative position at HH18 to the relative position at HH24 were found through linear regression of data obtained per arch (Fig. 5, Table S2). Variations within each stage are shown in the form of RSD (Fig. 5A). HH18 to HH24 change in peak WSS magnitude leads to significant arch artery correlations for PAAs IIIL, IIIR, IVR and IVL (Fig. 5B). The majority of these changes were associated with positive slopes, 30% were associated with a negative slope (Fig. 5B; Table S2). Because flow through these arches assumes a parabolic or Poiseuille flow profile, we can determine which factor (morphological or hemodynamic) drives such slope changes. As Poiseuille WSS magnitude is proportional to flow over radius cubed (Q/r³) or equivalently mean velocity over radius (v_mean/r), a positive slope suggests that change in flow or mean velocity are locally driving growth in the arches at these stages, whereas a negative slope implies the geometric feature (D_max, D_min, Area) is dominating the local growth between stages over flow or mean velocity. This analysis revealed that WSS played a significant role in promoting growth between stages, with WSS driving growth in IVL area changes, and PAA IIIL and IVR minimum diameter changes (Table S2). The same technique of quantifying a relationship between hemodynamic and geometrical changes was used for pressure in the form of pulsatility (difference between maximum and minimum pressure over a cardiac cycle). Pulsatility correlates with diameter changes in the early arch arteries (Fig. 5B). For HH18 to HH24 growth, significant trends were found for PAA IIIL. D_max was associated with a positive slope, whereas D_min was associated with a negative slope. Pulsatility can be considered a secondary driver of growth in HH18 to HH24 remodeling.

One of the strengths of the current technique is the ability to examine the changes in an individual (subject-specific configuration), and to highlight and learn from the changes in configurations that are at the same stage. Rather than assessing an overall general trend based on idealized composite geometries, this new technique allows for local as well as global variability analysis and therefore fosters an understanding of multiple forms of healthy development and how the progression comes about. Thus, as a first application of this technique, we examined the variations in HH26 embryos. Because of the major changes taking place in HH26 embryos, there is a larger bell curve of variation (Fig. S2), that convey a story of PAA growth and remodeling. These embryo to embryo variations in PAA geometry are especially pronounced during major transition phases where multiple configurations may be present at the same HH stage (Kowalski et al., 2013). Despite identical wing and limb bud characteristics, statistically significant differences were found for CSA and ellipticity parameters among subsets. HH26 marks a major time point of change, as cardiac septation or partitioning of the heart has begun. The aorticopulmonary septum has lengthened toward the heart and entered the distal OFT (Waldo et al., 1998), separating the OFT into an upper and lower section due to protruding mesenchyme right before the entrance to the PAAs. The large lengthening of the aorta and pulmonary trunk are reflected in Fig. 6C. Differences in timing and/or rate of lengthening, rotation and septation of the arches resulted in HH26 subsets that we refer to as ‘slow’ and ‘fast’ (Fig. 6A). The three configurations present at this stage allow for a better understanding of how the vessels evolve and the effect of OFT septation on the arches themselves. PAA cross-sections continue to become more round through HH26, obtaining an almost perfect circular cross-section at HH26-slow (Fig. 6A,B blue) before elongating in the vertical (cranial-caudal) direction (Fig. 6A,B red). HH26-fast embryos still feature this characteristic elongation in the vertical direction (Fig. 6B green). HH26-fast embryos exhibit more elliptical (particularly in the latter half) intricate arches with smaller CSAs when compared with the HH26 subset.

Owing to morphological changes in how arches connect to OFT (Fig. S4), flow is no longer required to mount the aortic sac in order to distribute to the more cranial arches as it was with the previous stages. The average flow distribution between arch pairs is 25% to 21% to 55% for PAA III to PAA IV to PAA VI arch pairs. The right and left sixth arches begin to receive the majority of flow, between 55% and 72% of the flow in four of the five HH26 embryos, and 25% of the flow in the less developed HH26-slow embryo. Pressure continues to increase substantially from HH24 to HH26. Across the HH26 embryos, peak pressure varies between 9350 dyne cm⁻² (7 mmHg) and 15,000 dyne cm⁻² (11.25 mmHg) at peak inlet flow (Fig. 7A), double that of the HH24 embryos. Despite the large increase in overall pressure magnitude, average pressure per arch, 6964 dyne cm⁻² (±214), is up only 15% from HH24 embryos, highlighting both the substantial increase in PAA length and the power of the septated OFT inlet to dissipate pressure. The caudal most arch pair (here, arch VI R and L) no longer consistently bears the largest pressure magnitude per arch length. In general, pressure is dissipated almost equally along OFT junctions and dissipates in a uniform fashion along the bilaterally paired PAAs. Pressure is distributed least uniformly in the HH26-slow embryo (HH26-5, Fig. S1C), where the L caudal VI arch still receives the largest pressure load, as with its HH24 counterparts. WSS levels remain at the same magnitude as those of HH18 and HH24 embryos, with overall WSS per arch magnitude of 91.48 (±21.56) dyne cm⁻², a slight decrease from HH24 embryos. Given these CFD obtained pressure, WSS and flow values, we again sought to determine whether WSS continued to drive these changes and whether pulsatility was a secondary driver. We performed linear regressions on HH24 to HH26 change in hemodynamic and change in morphological parameter data. As cardiac septation is a major programmed event, we hypothesized that the morphological changes would lead to the necessary hemodynamic changes. For peak WSS, significant correlations were found for PAA IIIR, IVR, IVL, VIL and VIR. As hypothesized, slopes for all HH24 to HH26 growth were negative, indicating that the geometric feature is dominating the local growth between stages over flow or mean velocity (Fig. 7B). Pulsatility maintains its role as a secondary driver of PAA growth with significant trends found for PAA IVL and VIL. Through the HH26 subsets, we see how the ‘lengthening and rotation of the aorticopulmonary septum toward the heart’ translates to a simultaneous pulling of the arch arteries from their OFT inlet and stretching along the cranial dorsal axis. These movements render the PAA cross-sectional areas more circular and then oval in the cranial-caudal direction, and work to maintain shear stress levels, while distributing the mounting pressure in a more uniform fashion. The effect of ellipticity on hemodynamic parameters is further detailed in Fig. S3.

To validate hemodynamic outputs from the computational simulations, in silico flow measurements and flow distribution of the PAA network were compared with experimental measurements (Fig. 8). At the level of the DoA, HH18 and HH24 flow curves are compared with those of Yoshigi et al. (2000), while HH26 average flow is compared with that of Hu and Clark (1989) (straight lines). Ultrasound velocities (Fig. 10) were converted into flow curves based on average PAA area along the arch and on Poiseuille or plug flow assumptions for the spatial profile. In general, there is a good agreement between computational and experimentally obtained values, with a significant difference between in vivo- and in silico-derived arch flow splits existing only for arch IIR based off a two-tailed Mann–Whitney U-test. When comparing in vivo and in silico flow splits across all arches, no significant changes are found for any of the present stages (two-tailed Mann–Whitney U-test). Although flow distribution hierarchy in terms of which arch received the most and the least amount of flow was consistently conserved between experimental and computational measurements, HH18 IIR and IIL flow was more constrained in silico than in vivo. Nevertheless, these comparisons indicate that in silico approximations of in vivo conditions are reasonable. Pressure waveforms were not measured or used as a direct input for boundary conditions in 3D simulations; rather, these values were indirectly controlled through the 0D bounds. Values are consistent with literature values measured at the level of the DoA by Yoshigi et al. (2000) and Zahka et al. (1989).

In an effort to understand the physical transformation occurring between stages, rather than change in stage to stage end points, statistical shape modeling was used to transform HH18 (using only the persistent arches, and not the regressing PAA II pair) to HH24 (without the growing PAA VI pair), and then HH24 (all 3 arch pairs) to HH26 (Fig. 9A). Statistical shape modeling (Deformetrica; Durrleman et al., 2014) was chosen as it handles complex amorphous structures and creates an unbiased indication of the most likely path of morphogenesis if overall distance minimization between shape distributions was the controlling variable (see Materials and Methods). Although the transformed geometries are not identical to the mean average geometry (template) at the transformed stage, understanding what displacement fields are necessary for embryos to approach the template shape leads to a better understanding of the changes taking place. HH18 to HH24 surface area displacement is highest in the cranial DoA branches and the ventral half of PAA IV. These patterns differ greatly from both the pressure and WSS HH24 maps (Fig. 4B). The magnitude of HH24 to HH26 3D surface displacement is more than double that of HH18 to HH24 deformation, consistent with the doubling of both the change in circulating blood volume (Fig. 7, top) and the change in pressure magnitude from HH24 to HH26. The HH26 displacement field resembles that of HH26 pressure maps (Fig. 7A), with displacement highest at the OFT inlet and dissipating through the arches toward the DoA. High displacement values also appear in the right cranial arch and the caudal-most section of the DoA on the transformed HH26 embryo. Although pressure (in the form of pulsatility) only played a small role in local diameter and area changes, pressure magnitude may play a larger role in overall arch movement and elongation from HH24 to HH26.

The second numerical model, which breaks down the PAA morphology into electrical analogs, offers another concrete way of examining growth by quantifying ease of flow transport through a vessel. From a flow perspective, each vessel or flow conduit, including the aortic sac, can be summarized through a single resistance value. The higher the resistance, the harder it is for flow to travel through and the more force is needed to ensure continuous circulation. To summarize how the system as a whole changed, resistance values for PAA manifolds were combined into a single resistance value per embryo, called the equivalent resistance (R_eq). There is a 7% increase between R_eq values for HH18 and HH24 embryos. A 106% increase exists between HH24 and HH26 R_eq values. Coupled with the displacement field metrics (Fig. 9A), the R_eq increases underscore the effects of vessel morphology on the functionality of the PAA system. On an individual arch artery resistance level, the general trend is a decrease in resistance values between HH18 and HH24 followed by a general increase in resistance values between HH24 and HH26. Statistically relevant changes exist only between HH18 and HH24 arch IVLs (P=0.035) and between HH24 and HH26 arch VIRs (P=0.045). When displayed in the context of a resistance versus area graph (Fig. 9C), regression, growth and maturation trends are consistent with generally accepted ideas of vascular growth, particularly as it relates to length and CSA changes (Fig. 9). The regressing PAA IIR and PAA IIL have the highest resistance and smallest mean CSA. HH18 PAA IVL and HH24 PAA VIL have almost identical resistance area values. These caudal arches are both growing in on their respective stages. For the established arches (PAA IIIR and PAA IIIL), HH24 resistance values increase from HH18 to HH24, consistent with the narrowing of the vessel diameters and decrease in CSA. The decrease in resistance value for maturing vessels is consistent with Poiseuille estimated resistance (1/Area² or equivalently 1/Radius⁴). Similarly, elongation that corresponds to an increase in length and decrease in area is associated with an increase in resistance. Resistance-area arch trends across stages (Fig. 9C) are fitted by an exponential function (R²=0.66). WSS versus surface area arch-specific values (Fig. 9C) resemble that of a third order polynomial (R² value=0.56), and indicate size-specific sensitivity zones for WSS-driven surface area changes. WSS-surface area graphs again underscore WSS maintenance of a 20-160 dyne cm⁻² range. Percent flow versus volume change graphs (Fig. 9C) highlight the evening out of the HH26 cranial most arches (flow percentage values, y-axis) and extreme increase in robustness/volume of the caudal ‘trunks’ (volume change values, x-axis). Percent flow-volume change graphs also emphasize how hemodynamic and morphological characteristics from one stage promote changes in the subsequent stage. These trends are consistent with information obtained from hemodynamic-morphology graphs (Figs 5 and 7). If hemodynamics were found to lead growth of a particular PAA from one stage to another, a small percentage flow value will lead to a small volume change (e.g. HH18 PAA IVR, IVL). If morphology were found to lead growth, a small percentage flow value may be associated with a large volume change (e.g. HH24 PAA VIR, VIL).

In order to determine which areas along the arch were the most tightly controlled in terms of hemodynamic and morphological parameters, we examined relative variability of pressure and WSS in conjunction with relative variability of D_min, Area and D_max. This local tracking is important, as regional differences in mechanical force, rather than overall levels, promote vascular remodeling. As morphological changes sometime lead hemodynamic changes, and vice versa, highlighting combined least-variability areas may help identify regional promotors of change. Deviation from mean was calculated for all arch sections among all embryos (ten sections per arch, five embryos per stage) (Fig. S5). Mean deviation was taken to be the difference between the mean and specific value normalized to the mean. The least variable region, from a combined hemodynamic-geometric perspective, was determined for each category (P-D_min, P-Area, P-D_max, WSS-D_min, WSS-Area, WSS-D_max) and each stage (HH18, HH24, HH26) (Table 1). The least variable regions for HH18 embryos were sections 1 and 2 along the arch, whereas the least variable regions for HH24 and HH26 embryos were in the latter half of the arch, particularly section 10. Surprisingly, trends did not always correlate with significant WSS- and P-driven area diameter changes, as determined through linear regression (Figs 4B, 6B, ST2). For example, WSS was found to be a major driver of HH18-HH24 Area IVL trends (P<0.0001, R²=0.86), but the least variable WSS-Area regions are found along PAA IIIL for both HH18 and HH24 embryos. Conversely, WSS was found to be a major driver of HH18 to HH24 PAA IIIL D_min changes (P=0.001, R²=0.76), and the least variable HH18 and HH24 WSS-D_min regions are found along this arch, as well as the least variable HH24 P-D_min region. We believe the highlighted regions should be closely watched when determining whether development is proceeding in a healthy manner. These regions should also be examined at the cellular/molecular level to determine whether their composition is the same as the surrounding vessels; they may serve as markers of normal development.

The dissemination of comprehensive three-dimensional atlases of development can serve as a powerful reference (de Bakker et al., 2016; Rana et al., 2014). More powerful still is the establishment of a range of healthy in addition to general trends for 3D morphogenesis. Examining variability to understand spatiotemporal morphogenetic and hemodynamic tolerance ranges, can help identify growth/adaptation mechanisms, determine risk corridors for malformation and predict locations of primary causation. RSD templates, like the one shown in Fig. 5A, help to identify which areas along the arch are more conserved and where high zones of variation (‘hotspots’, orange/red) exist. In HH18 embryos, the positions along the arch with the largest variation exists both along the IVL arch, in keeping with the fact that this vessel is growing inwards, and the regressing PAA IIL. With the growing in of PAA IVR, the highest relative variation exists along the center segment of the arch, suggesting an outward growth from the connections at either end. The beginning and end arch segments are also the most tightly controlled on a geometric-hemodynamic variability level (Table 1). As PAA IV (right and left) grows in, the vessel is narrowest through its midpoint section (Lindsey et al., 2015). The same is seen for PAA VI. As these vessels mature, the center segments obtain comparable area-diameter values with that of the rest of the vessel. Vessel regression follows an opposite trend, with apoptosis and reduction in vessel lumen diameter beginning at the midpoint before the proximal and distal ends disappear (Molin et al., 2002). HH24 embryos show the largest variation on right side vessel midpoints, as well as the first half of PAA IIIL. From HH18 to HH24, flow is shifting from flowing primarily through arch pair III to a more even flow distribution. The largest relative area variation for HH26 embryos takes place at the PAA IVR midpoint. By this time, the vessels have been separated into an ascending aorta and pulmonary trunk, with the most caudal VI arch pair receiving the majority of the flow and PAA IVR receiving the least amount of flow. Section 1 of PAA IVR has the least variable HH26 WSS-D_max zone, consistent with systematic flow percentage reduction to this arch.

WSS maintained a central role in PAA growth. The relatively constant magnitude of WSS across stages (20-160 dyne cm⁻²) suggests that this is a value to be maintained, in keeping with known trends (Beloussov, 2008; Taber, 2009). Linear regressions revealed a total of ten trends for HH18 to HH24 growth, and 13 trends for HH24 to HH26 growth. Twelve of these trends regulated growth on the right side while 11 of these regulated growth on the left side (Table S2, top). In addition to trends within arches, when reduced to a single value per arch, HH18 to HH26 WSS-surface area graphs follow a nonlinear function, monotonically increasing across stages (Fig. 9C). As hypothesized, a combination of pressure and WSS maps could account for many of the structural changes seen between stages over this critical window of development. Optimal WSS values may vary according to transmural pressure as determined by Pries et al. (1995). HH18 WSS zone of 80-160 dyne cm⁻², may therefore induce large area/diameter changes in HH18 embryos (mean pressure per arch 3387 dyne cm⁻²), but only slight deviations in HH24 embryos (mean pressure per arch 6964 dyne cm⁻²). Pressure plays a key role in the cardiovascular system, as it must be large enough to push the necessary amount of blood through the arches and downstream vessels in order to meet the embryo's metabolic needs (Kowalski et al., 2012; Pries et al., 1998; Zamir, 1977). In addition to the embryo's metabolic need for function and growth, the transport system itself has a metabolic cost (wall tunica media generation and blood volume maintenance) (Murray, 1926). PAA vessel morphology therefore results from an optimization process in which mechanical and metabolic stimuli trigger constant adaption (Kowalski et al., 2012; Murray, 1926; Pries et al., 1998; Taber, 1998).

To further characterize the effect/consequence of stage-specific variation in a 3D growth context, statistical shape modeling was used to transform the HH24 mean template into the HH26-slow and HH26-fast configurations. The overall displacement field magnitude remained the same for the HH24 to HH26-slow transformation at the level of the outflow tract (Fig. S6), which grew more in the vertical direction than the horizontal direction (unlike the longer more finessed HH26 and HH26-fast embryos). Although the HH26-slow outflow tract experiences the same magnitude of surface area displacement as the HH26 mean template embryo, the section of the arches most proximal to the heart experienced ∼25% less surface area displacement. The mean displacement across the entire HH24 template to HH26-fast template displacement increased by ∼21% (to 1.44). This increase was experienced more in the OFT than in the arches themselves. PAAs experienced a mean arch displacement of 0.73 (±0.16) cm between the large developmental transformation that is HH24 to HH26 growth. Their displacement maps support the theory of differential growth of the distal segments relative to the proximal.

A quantitative understanding of PAA growth and morphogenesis is necessary to uncover mechanisms of clinically significant outflow tract and great vessel malformations. These malformations involve substantial differential growth, rotations and motions of tissues formed from multiple cell lineages and in continuous engagement with complex hemodynamic environments. Purely genetic approaches lack the precision necessary, and even so quantitative tools for analyzing structure-function relationships across a population are sorely lacking. Here, we establish an integrated set of quantitative population-based morphological and hemodynamic analysis and simulation tools, and apply them to elucidate specimen-specific and group aggregate variations in flow, WSS and pressure. We quantified stage-specific vessel-specific functionality through resistance values. Results confirmed WSS as a major driver of PAA growth, and highlighted pressure as a secondary driver. Flow and pressure magnitudes increased dramatically across stages, while WSS magnitude was maintained. Side-specific flow dominance was not definite until HH26. Each stage was marked by a characteristic vessel shape that reflects the movement of the OFT and descending of the heart. By examining stage-specific characteristics as well as growth between stages from a number of perspectives, we define and quantify healthy PAA morphogenesis at this crucial stage of cardiac development. Ultimately, quantitatively interrogating variation in morphological outcomes is the best way to achieve new mechanistic insights into morphogenesis.

Fertile white Leghorn chicken eggs (Cornell Poultry Farm) were incubated blunt-side up for 3, 4 and 5 days in a continuous rocking incubator at 37.5°C. Embryos were randomly selected for each day endpoint for measurements and staging, geometric scanning, and modeling. This was intentional so as to study natural variation (that stems from genetic and epigenetic variability) within well-controlled conditions. In order to prepare embryos for geometric scanning, embryos were removed from the incubator at the appropriate stage and subsequently dissected away from their yolk sac. Injection micro-needles were fashioned from pulled capillary tubes (0.75 mm ID) cut to 20-35 µm inner diameter via a microforge (Glassworx). A micromanipulator (model M3301L, World Precision Instruments) was used to position the needle into the apex of heart. The embryo's vascular system was flushed with phosphate-buffered solution followed by 4% paraformaldehyde to preserve inner vascular volumetric integrity. The embryos were then left in 4% paraformaldehyde for 24-48 h before being transferred to a 70% ethanol solution and stored (possibly for several weeks). Embryos were brought up to a 30% ethanol solution before being transferred to a diluted form of Lugol solution, aqueous potassium iodide and iodine (Sigma-Aldrich, L6146). The embryos were soaked in Lugol's and the solution was changed over several days until the embryos no longer took up any iodine. Embryos were then dehydrated down to 100% ethanol, placed in polymerase chain reaction tubes and sent to undergo 3-4 µm nano-computed tomography scans (nano-CT).

Outflow tract (OFT) velocity and that of the three paired pharyngeal arch arteries were measured using B-mode guided (Movie 1, Fig. 10) Doppler Ultrasound (Vevo770 and Vevo 2100, Visualsonics). Warm Earle's balanced salt solution was used as an aqueous conduit between the embryo and the ultrasound scanhead. Embryos were kept at 37.5°C during imaging by being placed in a heated water bath as previously noted (Yalcin et al., 2010). After imaging, the embryos were transferred back to the incubator and allowed to recover. For each analyzed stage, 11-15 Doppler recordings were averaged to obtain a summary velocity profile for each vessel of interest.

0D electric analog representations of the arch artery network were created for HH18, HH24 and HH26 geometries. Sundials initial value problem solver Implicit Differential-Algebraic solver, IDA (Lawerence Livermore National Laboratory, Livermore) was used to obtain solutions to the 0D circuit. The full 0D model consisted of the arch artery system, including the aortic sac inlet, three arch pairs, cranial and caudal sections of the DoA, and the distal components associated with these outlets, according to the geometry resulting from the ‘In silico geometry preparation’ section. The distal components (Fig. S4, dashed boxes), an RCR or Windkessel model, were adapted from Yoshigi et al.'s single compartment 0D model of HH18 and HH24 embryo circulation, as seen from the PAAs (Yoshigi et al., 2000). The values of the arch arteries and flow conduits that blood traveled through to reach the outlets were obtained from the 3D simulation. Iterations between the 0D and 3D-0D simulations were necessary to tune the 0D boundary conditions (detailed below in ‘Multiscale 3D-0D simulations’), as shown by Pant et al. (2014). 0D representation of the HH18 and HH24 geometries are shown in Fig. S4. A separate 0D model exists for HH26 morphology where the onset of OFT septation requires that the aortic sac resistance be divided into two. Resulting pressures after this division are labeled Ps1 and Ps2. 0D abstraction of the 3D domain is represented in black; the Windkessel boundary conditions are represented in red and outlined with dashed lines (Fig. 2, Fig. S4). Electric circuit models were updated after 3D simulations according to Ohms law.

Embryo-specific 3D geometries of the HH18 (day 3), HH24 (day 4) and HH26 (day 5) PAAs were generated by importing nano-CT images into MIMICS (Materialise) and 3MATICs (Materialise). PAA geometries extended from the distal outflow tract to the DoA and paired cranial aortae (Fig. 2). Geomagics Studio 10 (Geomagic) was also used for the preparation of 3D geometries for CFD. All embryos were scaled to account for the difference between dehydrated and native vessel size, as determined by a series of in vivo diameter verifications. India ink and Texas Red Dextran were used to obtain native vessel size across stages, results were compared with that of 3D reconstructions and a scaling factor was generated.

Blood was treated as a Newtonian fluid with constant hemodynamic properties (ρ=1060 kg/m³, μ=3.71×10⁻³ Pa.s) and rigid, impermeable vessel walls were assumed with no slip boundary conditions. Flow was simulated using the in-house code FELiScE (finite elements for life sciences and engineering; gforge.inria.fr/projects/felisce; Pant et al., 2014) on a high-resolution unstructured Cartesian grid with finite-element numerical treatment. For 3D mesh generation and adaptation, feflo and ghs3d (Loseille and Löhner, 2010) were used. Grid sensitivity analysis was conducted on a control PAA model for each stage in order to ensure consistency and reliability of the numerical solutions for all simulations presented in this study, beyond which resulting mass-flow redistributions were insensitive to further Cartesian grid refinements.

The Vascular Modeling Toolkit (vmtk), supported by Orobix Srl, was used to obtain cross-sectional area and diameter measurements taken along the centerline of each geometry (Piccinelli et al., 2009). Once the centerline was obtained, normal and tangent functions were calculated for one in three points along the centreline, in order to obtain the equation of a plane normal to the centerline at each point. Each section was imported into Ensight and its orientation verified before calculations were performed. Ensight (Computational Engineering International), was used for all post-processing, visualization and measurement of flow properties.

Pulsatility was taken to be pressure at peak flow minus the pressure at the minimum flow (P_max−P_min). Pressure drop was taken to be average pressure in the most proximal end of a segment minus average pressure in the distal end of a segment. An equivalent resistance value (R_eq) was calculated across stages according to Eqn 1. P_in was taken to be pressure at the OFT. P_out was taken as a combination of the three outlet pressures (left cranial, right cranial and caudal dorsa aorta outlets) multiplied by their respective flow distribution.

Deformetrica (Durrleman et al., 2014) statistical shape modeling software was used to create stage-specific average PAA surfaces and to grow one stage-specific cohort to the next. The averaged geometrical surface served as a template for mapping morphological characteristics to the surface for each of the HH18, HH24 and HH26 geometries. Registration was also performed between the HH18 and HH24 templates, as well as between the HH24 and HH26 templates. The HH18 template was then deformed towards the HH24 template geometry and the HH24 template deformed towards the HH26 template through a combination of distance minimization and transformation algorithms driven by moment vectors over the ambient space. Through the statistical shape modeling software, geometrical surfaces do not need point-by-point mesh correspondence, and their dissimilarity is evaluated by a metric on a mathematical current (flux of a vector field) framework that deforms the underlying 3D space, rather than on anatomical landmarks themselves. This non-parametric technique allows for the handling of complex amorphous structures, and is particularly well suited for congenital heart defect models (Guibert et al., 2014). At the basis of this non-parametric comparison are currents. Currents characterize a shape as a distribution of shape features. Surfaces are embedded in 3D space and shape variations are measured based on deformations of the underlying 3D space. Shapes are compared by computing how distributions of features differ. In the ‘forward approach’, which is applied by this model, each subject is a deformation of the anatomical mean shape (the template) plus a residual (shape features not represented by the template). The deformations are described by diffeomorphisms or functions that preserve the topology of an object, mapping one manifold (the template) to another in a continuous fashion. Computational surface meshes serve as inputs into the model. A mean template is iteratively computed as an average of all currents. At each iteration, the present template is registered to each of the study's subjects, computing the individual deformation. A new template is computed through a deformation in order to minimize the distance between the (newly) deformed template and target shape object (individual subject). A transformation function maps the template towards each individual subject shape, allowing each individual subject to be characterized by a multitude of deformation vectors rather than its actual 3D surface. This technique is called large deformation diffeomorphic metric mapping.

Morphological and hemodynamic changes were compared qualitatively and quantified when possible. Results were summarized in the form of mean and s.d. values over the course of one cardiac cycle. Paired t-tests were used where appropriate with P<0.05 denoting significance. Linear regressions and exponential decays were performed.

Statistical comparisons were made using GraphPad Prism. Linear regressions were plotted on hemodynamic versus geometrical parameter graphs for each of the HH18, HH24 and HH26 stages. The Friedman test was used to test for significant differences between HH26 arch subgroups. The Mann–Whitney U-test was used for flow split validation. Friedman's test is a nonparametric test used to compare three or more paired groups and by ranking the values within groups and then providing a one-way analysis of variance.

In order to summarize local stage-specific features, local RSD and deviation from mean values were calculated based on a point-by-point comparison throughout the arches. RSD was taken to the specific value divided by the average value across embryos (specific value/mean value). Deviation from the mean was taken to be the specific value minus the mean normalized to the mean (specific value−mean value)/mean value. Values for each arch were taken at ten locations along the arch, from beginning to end. In this way, values were compared at relative position, rather than at a precise distance for both morphological and hemodynamic parameters.

We thank the Cornell BRC Imaging facility, Pascal de Santa-Barbara (critical reading of manuscript), Jessica Ryvlin (ultrasound) and Alex Chang (ultrasound).