Vortex Dynamics in Trabeculated Embryonic Ventricles

Proper heart morphogenesis requires a delicate balance between hemodynamic forces, myocardial activity, morphogen gradients, and epigenetic signaling, all of which are coupled with genetic regulatory networks. Recently both in vivo and in silico studies have tried to better understand hemodynamics at varying stages of veretebrate cardiogenesis. In particular, the intracardial hemodynamics during the onset of trabeculation is notably complex—the inertial and viscous fluid forces are approximately equal at this stage and small perturbations in morphology, scale, and steadiness of the flow can lead to significant changes in bulk flow structures, shear stress distributions, and chemical morphogen gradients. The immersed boundary method was used to numerically simulate fluid flow through simplified two-dimensional and stationary trabeculated ventricles of 72, 80, and 120 h post fertilization wild type zebrafish embryos and ErbB2-inhibited embryos at seven days post fertilization. A 2D idealized trabeculated ventricular model was also used to map the bifurcations in flow structure that occur as a result of the unsteadiness of flow, trabeculae height, and fluid scale (Re). Vortex formation occurred in intertrabecular regions for biologically relevant parameter spaces, wherein flow velocities increased. This indicates that trabecular morphology may alter intracardial flow patterns and hence ventricular shear stresses and morphogen gradients. A potential implication of this work is that the onset of vortical (disturbed) flows can upregulate Notch1 expression in endothelial cells in vivo and hence impacts chamber morphogenesis, valvulogenesis, and the formation of the trabeculae themselves. Our results also highlight the sensitivity of cardiac flow patterns to changes in morphology and blood rheology, motivating efforts to obtain spatially and temporally resolved chamber geometries and kinematics as well as the careful measurement of the embryonic blood rheology. The results also suggest that there may be significant changes in shear signalling due to morphological and mechanical variation across individuals and species.


Introduction
Since the developing heart and embryo continue to function for some time in the absence of erythrocytes, it appears that the function of the early embryonic heart is not for the purpose of nutrient transport [1,2]. Instead, recent work suggests that the heart's function is to aid in its own growth [1,[3][4][5]. Two important roles for intracardial fluid dynamics in terms of proper cardiogenesis are to exert hemodynamic forces onto the ventricular lining and to advect morphogens [2,6]. These two fluid effects help regulate and drive organogenesis in developing embryos. Both shear stress and pressure may be key components to activating developmental regulatory networks by acting on cardiac cells [7] through a process called mechanotransduction. In this case, mechanical stimuli are transmitted by intracellular signalling pathways to the interior of the cell. Moreover, increased receptor-ligand bond formation may appear near the endothelial lining in regions of higher vorticity [8], which gives rise to greater mixing of chemical morphogens. These chemicals may act as epigenetic signals, which are then advected throughout the chamber [9,10]. It is clear that irregular hemodynamics leads to cardiomyopathies or embryonic death [3,[11][12][13][14].
The first heart beat occurs around ∼1 dpf when its morphology resembles a simple valveless tube. It is composed of both an outer myocardial and endocardial layer of cells. Around 1.5 dpf cardiac looping and chamber ballooning begins, where the heart drastically begins to reshape itself into a multi-chambered pumping system. Two chambers can be distinctly seen by 2 dpf, while endocardial cushions, the precursors to valve leaflets, are in the process of forming. Trabeculae, irregular muscular protrusions that expand from the inner surface of the ventricle, begin to form around 3 dpf. These stages are illustrated in Figure 1. Illustrations of significant stages of zebrafish heart development. At 1 dpf the heart begins pumping as a valveless heart tube, at 1.5 dpf cardiac looping and chamber ballooning begins, at 2 dpf the heart is a growing two-chambered pumping system with endocardial cushions forming, and at 3 dpf ventricular trabeculation occurs. The blue arrows indicate blood flow while the black dotted arrows represent mechanical stresses re-configuring the heart. Prior to trabeculation, the endocardial ventricular cells are smooth and polygonal in shape. During the onset of trabeculation, several endocardial cells become elongated and some extend cellular projections. Moreover, these cells appear slightly more depressed than the surrounding endocardial cells. The depressions progressively become deeper and wider such that the endocardial cells invaginate the cardiac jelly and extend toward the basal surface of the myocardium. Eventually, myocardial cells separate due to the potent endocardial cell invasion, and definitive trabeculae are formed [15]. Hence, trabeculae are composed of both myocardial and endocardium components.
Proper trabeculation requires well-coordinated cardiac contraction [16] and is particularly sensitive to local changes in the fluid environment [17]. It is thought that the trabeculae may serve as mechanotranductive structures and alter intracardial flows in a way that regulates shear stress and mixing near the endocardium [18][19][20]. Even if subtle trabeculation irregularities were masked, cardiac defects would magnify over time because of their effect on morphogenetic processes. For example, zebrafish embryos designed to lack normal trabeculation (ErbB2-inhibited) displayed severe cardiovascular defects including bradycardia (decreased heart rate), decreased fractional shortening, and impaired cardiac conduction [21]. Lack of trabeculae or irregularly formed trabeculae will cause irregular patterns of shear stresses. This in turn can cause dysfunctional myocardial activation patterns that are known to cause arrhythmias, abnormal fractional shortening, and even ventricullar fibrillation [22].
During the onset of trabeculation, the underlying fluid dynamics are particularly interesting due to the balance of inertial and viscous forces. The Reynolds number, Re, is a dimensionless number that describes the ratio of inertial to viscous forces in the fluid. It is given by Re = (ρUL)/µ, where µ and ρ are the dynamic viscosity and density of embryonic blood, respectively, and L and U are characteristic length and velocity scales. The characteristic velocity is often chosen as the average or peak flow rate, while L is often selected as the diameter of the chamber or vessel. When trabeculation begins as cardiac looping and ballooning progress, the Re is approximately 1. At this fluid scale, a number of important fluid dynamic transitions can occur. One notable feature is the transition to vortical (disturbed) flow and hence changes in flow direction. This transition is sensitive to the growing complex morphology, effective viscosity of the blood, and unsteadiness of the flow.
Disturbed blood flow patterns have been observed during heart development [23]. These flow patterns generally induce lower wall shear stresses (WSS) than smoothly streaming (non-disturbed) blood flow. Both contribute to remodelling and endothelial cell activation in different ways [24,25]. It has been shown that disturbed (vortical) flow patterns upregulate the expression of certain genes, such as Notch1, in endothelial cells during development [26]. As the heart undergoes dramatic morphological transformations, transitions to vortical and disturbed flow patterns may help guide morphogenesis through changing patterns in WSS or through other mechanotransductive mechanisms such as flow sensing through primary cilia [27][28][29]. Note that intracardial flows are both temporally and spatially varying such that the distribution of WSS is not uniform along the endothelium [19,20]. Hence, mechanotransducers will exhibit different responses, leading to differentiated cellular behavior [30]. Furthermore, as the heart grows, blood flow also increases [31]. The formation of complex structures along the ventricle, like trabeculae, may provide regions where disturbed flow develops which could lead to higher kinetic energy dissipation. This energy dissipation may facilitate proper ventricle contractile function and trabecular organization [20].
Due to the complexity of the cardiogenesis and the challenges of measuring flow patterns precisely, computational fluid dynamics (CFD) has become a premier tool for resolving the flow in embryonic hearts [19,20,[32][33][34][35][36][37][38][39][40][41]. For example, Liu et al. [34] simulated flow through a three-dimensional model of a chick embryonic heart during stage HH21 (after about 3.5 days of incubation) at a maximum Re of about 6.9. They found that vortices formed during the ejection phase near the inner curvature of the outflow tract. In 2013, Lee et al. [38] performed 2D simulations of the developing zebrafish heart with moving cardiac walls. They found that unsteady vortices develop during atrial relaxation at 20-30 hpf and in both the atrium and ventricle at 110-120 hpf.
More recently, Vedula et al. [19] and Lee et al. [20] used light-sheet fluorescent microscopy and reconstructed a 4D moving ventricle on which they based their CFD model. They were able to quantify spatially-and temporally-varying WSS along trabceular ridges (trabeculae "heads") and groves ("intertrabecular regions"). In particular, Lee et al. discovered that pulsatile shear-stresses developed along the ridges at 3 dpf in wildtype (WT) zebrafish embryos, while oscillatory shear-stresses (OSS) developed in the groves around 4 dpf [20]. Around 4 dpf, vortical flow patterns may be present within the intertrabecular spaces. Moreover, OSS were found to be substantially less at the trabecular heads, suggesting that OSS may be a possible regulatory control during cardiogenesis [19,20]. They also investigated differences in WSS between wildtype and mutant zebrafish hearts. The mutants they considered were ErbB2-inhibited zebrafish (suppresses trabeculation), gata1a morpholinos (lowers blood viscosity), and wea mutants (lower cardiac contractility). They found that total WSS was comparable in the WT and ErbB2-inhibited zebrafish; however, the gata1a morpholinos and wea mutants expressed significantly less total WSS. Another study, Battista et al. [40], found that trabeculae morphology has a significant effect on intertrabecular vortex formation, as does the presence of hematocrit and fluid scale. However, their study did not include an analysis of WSS, but instead referred to the tangential, normal, and total force magnitudes as potential proxies for WSS, although they observed similar trends to the spatially-averaged WSS over the course of a heart cycle.
The numerical work described above and in vivo measurements of blood flow in embryonic hearts [3,42] supports that vortex formation is sensitive to changes in Re, morphology, and unsteadiness of the flow. Santhanakrishnan et al. [35] used a combination of CFD and flow visualization in dynamically scaled physical models to describe the fluid dynamic transitions that occur as the chambers balloon, the endocardial cushions grow, and the overall scale of the heart increases. They found that the formation of intracardial vortices depends upon the height of the endocardial cushions, the depth of the chambers, and the Re. Their study only considered steady flows in an idealized two-dimensional chamber geometry with smooth, stationary walls.
In this paper, we present complementary studies to both Santhanakrishnan et al. [35] and Lee et al. [20] with the goal of revealing the bifurcations in flow structures that occur as a result of the unsteadiness of the flow, trabeculae height, and Re. First, we investigate the differences in the cardiac fluid dynamics between WT and ErbB2-inhibited (namely ErbB2 st61 and ErbB2 st50 ) mutants, to explore how vortex formation in the intertrabecular regions is sensitive to differences in morphology. We quantify the intertrabecular flow patterns mentioned (but not shown) in Lee et al. [20]. Next, we use an idealized geometry, based upon that of Santhanakrishnan et al. [35], to systematically sweep a parameter space consisting of trabeculae size, fluid scale, and unsteady flow effects to quantify fluid dynamics transitions.
It is important to note that, in this study, we perform 2D simulations using stationary boundaries. This permits direct comparison to the steady flows moved through the simplified chambers considered by Santhanakrishnan et al. [35]. It also permits exploration of a wide parameter space that covers the range of biological diversity, including other trabeculated vertebrate embryonic hearts and trabeculated invertebrate hearts. As such, the results highlight parameter spaces where bulk flow patterns are highly sensitive to small changes in morphology, effective viscosity, and unsteadiness of the fluid. In such parameter regions, it is critical to obtain highly resolved descriptions of morphology, pumping kinematics, and rheology of the blood. The simulation results shown here are not as realistic as those presented for the zebrafish embryonic hearts simulated in 3D with moving boundaries by Lee et al. [20] and Vedula et al. [19]. Our goal, however, was to map the parameter space sensitive to small changes in a relatively simple model. We believe that the results will serve as motivation for more detailed three-dimensional studies. We also argue that the stationary boundaries and 2D approximations serve as a reasonable starting point given the relatively low Womersley number, Wo ≈ 0.77, since unsteady effects become significant for Wo > 1.

Methods
Two-dimensional computational fluid dynamics (CFD) simulations were used to quantify the flow fields within a biologically-realistic and an idealized model of a trabeculated ventricle of the zebrafish embryonic heart. We will discuss the model geometry construction below and emphasize that the purpose of this study is to illustrate that flow patterns in and near trabeculated regions are particularly sensitive to small changes in morphology and scale in a biologically relevant parameter space. For a detailed discussion on the numerical method used to solve the fluid-structure problem, see Appendix A.

Embryonic Zebrafish Model Geometry
We estimated the structure of trabeculae from in-vivo image data taken from embryonic zebrafish for the purpose of simulating flow through realistic ventricle geometries. Five images of stained cross sections of both wildtype and mutant embryonic zebrafish were taken from Liu et al. [21]. The mutant and transgenic lines used were ErbB2 st61 and ErbB2 st50 . Figures from the original paper were cropped out and manually traced by recording pixel locations along the boundary with the open-source Python package Argus [43]. Smooth-spline approximations of the trabeculated ventricle walls were generated and sampled to generate the discretized computational geometry as seen in Figure 2. Figure 2 shows the computational geometries extracted from microscopy images [21]. Note that the images were acquired using a Nikon Te-2000u microscope (UNC-CH, Chapel Hill, NC, USA) at a rate of 250 frames per second with a high-speed CMOS camera (MiCam Ultima, SciMedia, UNC-CH). The geometries chosen for this study were of wild-type embryos at 3 and 5 dpf and an ErbB2 inhibited embryo at 7 dpf. Note that, at 5 dpf, ErbB2-inhibited embryos show little to no signs of trabeculation [21]; however, at 7 dpf, there is the onset of trabeculation, as shown above. Furthermore, we considered other wild-type embryos at 80 hpf and 5 dpf to span biological variation; see Figure A1 in Appendix B for their geometries. Note that we distinguish between both 5 dpf embryos by giving them an (I) or (II) designation.  [21]. There are two wild-type (WT) zebrafish cases shown (3 dpf and 5 dpf (I)) as well as the geometry taken from an ErbB2-inhibited zebrafish at 7 dpf.
To ensure the simulations were conducted at the appropriate fluid scale, we computed the biologically relevant Re to encapsulate the correct fluid dynamic regime. For the Re within the ventricle of a 4 dpf wild type zebrafish, the characteristic velocity, V z f , was taken as the average of the minimum and maximum velocity measured in vivo, and the characteristic length, L z f , was taken as the width of the ventricle from Figure 3a. The Re was then calculated as where V z f = 0.75 cm/s [3], ρ z f = 1025 kg/m 3 [5], µ z f = 0.0015 kg/(ms) [44], and L z f = 208 µm.
The characteristic frequency was chosen as f = 3.95 beats/s [18]. The dimensionless frequency may then be calculated asf The simulations using the biologically realistic geometry were performed at Re ∼ 1. The length scale was taken directly from the images of the embryos and the velocity, frequency, and kinematic viscosity (the ratio of the fluid's dynamic viscosity to density) were taken from the literature as described above. Figure 4 shows how the geometric parameters were measured from the reconstructed geometries. Note that, although there is geometric variation between each case, the parameters are labeled consistently. We do not, however, vary the Re for the simulations that use the realistic geometry. Note that we do vary Re for the idealized trabeculated ventricle geometry described in Section 2.2.  [21]. This snapshot was taken immediately before the systolic phase. The protrusions into the ventricle are trabeculae. The image was taken from Tg(cmlc2:dsRed)s879; Tg(flk1:mcherry)s843 embryos expressing fluorescent proteins that label the myocardium and endocardium, respectively [21]; (b) simplified diagram showing the basic idea behind our idealized geometry. Blood flows from the atrio-ventricular canal into the ventricle and then proceeds into the bulbus arteriosus; (c) further idealization of the computational model which is now flattened. The geometric parameters are as follows: a V and b V are the semi-major and semi-minor axis of the elliptical chamber, h T and r T are the height and radii of the trabeculae, and w AV and w SV are the widths of the atrioventricular (AV) canal and sinus venosus, respectively.

Idealized Model Geometry
A simplified two-dimensional geometry of a 4 dpf zebrafish's trabeculated ventricle was constructed using Figure 3a, from Liu et al. [21]. The ventricle was idealized as a half ellipse, with semi-major axis a V and semi-minor axis b V . It is tangentially laid within a channel, which emulates a cavity-flow geometry. The channel models the atrioventricular canal (AV canal), with width w AV set equal to the width of the sinus venosus (SV), w SV . Six equally-spaced trabeculae were aligned within the ventricle. The idealized model geometry is illustrated in Figure 3c. Note that the aspect ratio of the chamber was selected for direct comparison with Santhanakrishnan et al. [35].
The idealized trabeculae geometry was modeled using the following perturbed Gaussian-like function, where r T and h T are the radii and height of the trabeculae, respectively. The full idealized geometry can be seen in Figure 3c. The geometric model parameters are given in Table 1. All parameters were scaled from measurements taken from Figure 3a. The parameters describing the ventricle were held constant and are given as the chamber height, b V , chamber width, a V , and width of the AV canal and SV, w AV and w SV respectively. Note that the radii of the trabeculae, r T , was held constant in all numerical simulations, while the height of the trabeculae, h T , was varied.

Parameter
Value The temporal parameter values were chosen to keep the dimensionless frequency fixed at 0.10 for the pulsatile simulations, analogous to the model presented in Section 2.1. The Re was varied by changing the kinematic viscosity, ν = µ/ρ. The computational parameters are found in Table A1. For the simulations, the Re sim is calculated using a characteristic length of w AV and a characteristic velocity set to V in (steady inflow) or 1 2 V in (pulsatile inflow). The simulations were performed for Re sim = 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 20, 30, 40, 50, 100. The stiffness of the tether springs were chosen the minimize the deformations of the boundary to below 1% of the chamber diameter. This approach is used in the immersed boundary method to describe a nearly rigid boundary.
Note that we consider Re two orders of magnitude higher and lower than the relevant Re ∼ 1 for a 4 dpf zebrafish heart. The reason for this is twofold: (1) to map out the parameter space outside of the biologically relevant range and (2) to provide insight into fluid flows in other types of trabeculated hearts, such as those of some invertebrates [45,46].

Results
Below, we present the flow patterns and velocities for both biologically realistic and idealized 2D models of trabeculated ventricles. The cases with biologically realistic geometries were of WT zebrafish and an ErbB2-inhibited mutant to contrast the intracardial and intertrabecular fluid dynamics of an embryonic zebrafish heart during development. The idealized geometry was used to systematically sweep over a parameter space to describe transitions in flow patterns for pulsatile flows, changes in trabecular height, and Reynolds Number Re. In the idealized geometry case, the Re was varied from 0.01 to 100, and the trabecular heights were varied from zero to twice the biologically relevant height. We also quantified flow for both steady and pulsatile cases.
Streamlines were used to show the path that a passive particle would take in the flow. The streamline graphs were generated using the VisIt visualization software (v. 2.7.2, Lawrence Livermore National Laboratory, Livermore, CA, USA) [47]. The streamlines are drawn by making a contour map of the stream function, since the stream function is constant along the streamline. The stream function, ψ(x, t), in 2-D is defined by the following equations: The streamline colors correspond to smooth, streaming flow (blue) and vortical flow (orange). Figure 5 gives the flow field streamlines for the steady inflow cases with biologically realistic geometries at Re ∼ 1. Parabolic inflow was used that accelerates from rest to a constant velocity, as detailed in Appendix A. The simulation data is given once the flow has reached its steady state. Results are shown for a 3 dpf, 80 hpf, and two 5 dpf wild-type embryos as well as an ErbB2-inhibited embryo at 7 dpf. In all cases, vortex formation occurred in the intertrabecular regions along the side opposite that to the sinus venosus (SV) and was within the well pronounced intertrabecular grooves. No intracardial vortices were observed and the flow smoothly steams from the atrioventricular canal to the sinus venosus.  Figure 6 illustrates that, in the regions of significant vortex formation, e.g., in the intertrabecular regions, velocities are much lower than in the intracardial region of the ventricle. In all cases, the velocities on the side opposite the SV experience much faster velocity decay towards the ventricular lining, whereas regions opposite the AV canal experience slower decay. Note that, in the 3 dpf and 7 dpf ErbB2 inhibited cases, the flow velocities along the line drawn from the AV canal to the intertrabecular region that is closest to the SV decrease in magnitude before increasing again and finally decaying to zero near the cardiac wall. This is due to the fact that the velocity is measured close to the ventricular wall that extends from the AV canal. Flow velocities measured between the trabeculae are significantly lower than those measured in within the middle of the chamber. The magnitude of velocity (cm/s) is given by the colormap, and the arrows show the direction of flow once steady state is reached (left). The panels on the right show the magnitude of velocity along lines drawn from the center of the AV canal to the ventricular lining between the trabeculae. Note that these lines are shown in the middle panels.

Pulsatile Flow through an Embryonic Zebrafish Heart
In this subsection, we show the results of the pulsatile inflow simulations where biologically realistic geometries are used at Re ∼ 1 and the frequencies are varied. Temporal snapshots of the streamlines are presented in Figure 7. The pulsatile inflow condition is detailed in Appendix A. Note that the period of the pulsation cycle is T. At the beginning and end of each pulsatile cycle, the flow is near zero. During the middle of the pulsation cycle, flow is maximal. Although this model is an idealization, we consider the first 50% T as analogous to diastole and 50-100% T as systole. All data shown is taken from the last pulse cycle such that periodic steady state has been reached. The geometries simulated include a 3 dpf, a 80 hpf, and two 5 dpf wild-type embryos as well as an ErbB2 inhibited embryo at 7 dpf. Vortices formed in all cases within the intertrabecular regions along the side opposite that to the sinus venosus, similar to the steady cases illustrated in Figure 5. The vortices are, however, dynamic; they change shape and size during a single pulsation cycle. Some additional vortices appear in the intertrabecular regions that were not present in the steady case. For example, consider the 80 hpf and ErbB2-inhibited cases where vortices appear on the side opposite the AV canal. Furthermore, large intracardial vortices appear between pulsation cycles in all cases when there is minimal inflow. Figures 8 and 9 give temporal snapshots for the 5 dpf WT (I) and the ErbB2-inhibited cases, respectively. The colormap illustrates the magnitude of velocity (cm/s) with accompanying velocity vector fields. The velocity field is consistent with Figure 7 and shows vortex formation. These results also illustrate the spatial gradient in velocity within the intracardial to intertrabecular regions. It is clear that, although the fastest flow moves towards the middle of the chamber, the velocity significantly tapers off by the time it reaches the ventricular lining, the region between trabeculae on the right-hand side of the ventricular geometries in Figures 8 and 9.   10-12 illustrate the magnitude of the velocity along lines drawn from the center of the AV canal to the intertrabecular regions for cases of a 3 and 5 dpf WT (I) and 7 dpf ErbB2-inhibited zebrafish geometries. Appendix C presents this data in a non-logarithmic scale, to complement the data given in Figure 6 in Section 3.1. Note that the 3 dpf WT and 7 dpf ErbB2-inhibited zebrafish exhibit smaller trabeculae than the 5 dpf WT (I) case.
As expected in all cases, peak flows occur within the chamber and decay as one moves toward the cardiac wall. Even in the steady inflow cases, the velocity decay is similar. The velocity decreases close to geometrically towards the intertrabecular region on the opposite side to the SV for the smaller trabeculae (3 dpf WT and 7 dpf ErbB2-inhibited). For the 3 dpf WT embryos in regions with pronounced trabeculation, the velocity decay does not monotonically decrease; it increases slightly within the trabeculae before decreasing to zero as one moves towards the endocardium. Similar non-monotonic behavior is seen for some intertrabecular regions farther from the SV in the 3 dpf WT and 7 dpf ErbB2-inhibited embryos. We also note that the magnitude of flow decreases much more rapidly for the 5 dpf WT (I) with pronounced trabeculation, indicating lower WSS within the intertrabecular regions.   Lastly for the case of pulsatile flow through realistic geometries, we computed the spatially-averaged wall shear stress (WSS) over the entire ventricular geometry, each pronounced trabeculae, and intertrabecular region, see Figure 13, as well as the oscillatory shear index (OSI). The OSI is a dimensionless metric of how aligned the WSS vector is with the time-averaged WSS (TAWSS) vector throughout one cardiac cycle. Recently, it has been suggested to be a possible regulatory mechanism of trabeculation itself [19]. Details regarding how WSS, TAWSS, and OSI were calculated can be found in Appendix D. Figure 14 gives spatially-averaged magnitude of WSS over the trabeculae, intertrabecular regions, and entire ventricle for the 3 dpf WT and 5 dpf WT (I) zebrafish geometries. The intertrabecular regions experience slightly less (∼10%) less spatially-averaged WSS than pronounced trabeculae. The majority of the WSS is felt elsewhere within the ventricle, e.g., near the AV-canal and outflow tract. See Appendix D.1 for a description of which trabeculae and intertrabecular regions were used for the analysis.   Figure 15 gives the OSI over the entire trabeculated ventricle for the 3 dpf WT, 5 dpf WT (I), and 7 dpf ErbB2-inhibited zebrafish geometries. High OSI ∼ 0.4 occur in both intertrabecular regions and on the trabeulae themselves. This is similar behavior to [19,20] who reported high OSI values in intertrabecular regions as well; however, in our simulations, not all intertrabecular regions nor trabeculae have high OSI. In particular, those nearing the outflow tract (top of geometry) experience less OSI versus those opposite the lower end of the AV canal.  Figure 16 shows the flow field streamlines for the case of steady flow through an idealized trabeculated embryonic ventricle. The inflow condition is detailed in Appendix A. The numerical simulations span five orders of magnitude of Re, varying from 0.01 to 100, while trabeculae heights were set to 0 ≤ h T b V ≤ 0.16. Note that the biologically relevant case is h T b V = 0.08 In the case of no trabeculae (left column), we find vortex formation only occurs for Re ≥ 15, in agreement with the findings of [35]. Moreover, this result appears consistent with fluid dynamics literature on transitions to vortical flow via a channel with an expanded region [48,49]. Shen et al. [48] and Mizushima et al. [49] both use rectangular cavities; however, similar transitions to vortical flow occur at consistent Re. For Re ≤ 10, the flow bends around the cavity and no flow separation occurs, that is, there are no transitions to vortical flow patterns. Hence, fluid mixing is not enhanced at for Re ≤ 10 and flow reversal does not occur. As Re is increased to 20, flow reversal occurs and a closed vortex is present along the left side of the cavity. The stagnation point is located between the orange and blue streamlines. To the left of this stagnation point, the flow moves along the endocardium from the right to left. To the right of the stagnation point, the flow moves right to left. As Re is further increased, the stagnation point moves to the right, and the intracardial vortex becomes larger until it becomes as large as the cavity itself for Re ∼ 100.

Steady Flow through Idealized Trabeculated Chambers
When half-size biologically relevant trabeculae are introduced into the model ( h T b V = 0.04), similar flow fields emerge for the case of Re ≤ 10. Although geometric perturbations now exist along the cavity lining, no flow separation occurs, whether intracardially or intertrabecularly, illustrating no vortical flow patterns nor enhanced fluid mixing at these trabeculae heights for Re ≤ 10. For Re = 20, we see a similar intracardial vortex to the case without trabeculae; however, this vortex weaves along regions with trabeculae. Furthermore, there is an emergence of an independent closed vortex along the right side between two trabeculae. For Re ≥ 50, we find the presence of one large intracardial vortex wrapping around each trabeculae.
For biologically relevant trabeculae heights, there are closed intertrabecular vortices for Re as low as 0.1, while no intracardial vortices are present at these lower Re. Shen et al. [48] saw a similar phenomenon with the formation of two vortices in the corners of their rectangular cavity. This is consistent with the formation of vortices near the only bottom of the trabeculae. On the other hand, interestingly, not all intertrabecular regions have closed vortices. As Re is further increased from Re = 5 to Re = 10, the intertrabecular vortices grow in size. As in previous cases, a larger intracardial vortex forms at Re = 20. On the left hand side of the cavity, there is smooth flow from left to right around the trabeculae. On the right side of the cavity, independent closed vortices form between the trabeculae, and the flow is from right to left. For Re ≥ 50, a large intracardial vortex forms and no intertrabecular vortices persist.  Furthermore, for the biologically relevant case of Re = 1, Figure 17 illustrates the magnitude of velocity from the intracardial center to the ventricular lining for various intertrabecular regions and trabeculae heights. It is clear that, for larger trabeculae heights, the velocity decays moving away from the intracardial center at a faster rate than those cases with shorter or no trabeculae.
When the trabeculae are at biologically relevant heights or larger, the velocity decays from the middle of the ventricle until a distance of that height away from the ventriclular lining. Within that trabecular valley, the velocity actually increases before further decaying when measured closer to the ventricular lining. This shows that a local minima in velocity magnitude exists, suggesting that trabeculae morphology plays an important role in governing the fluid dynamics nearing the ventricular lining. This supports that there is intertrabecular vortex formation in these cases, as shown in Figure 16. The velocities measured in such intertrabecular regions are less than those shown in the cases with no or smaller trabeculae. Figure 17. Measurements of the non-dimensional magnitude of velocity are provided for the Re = 1 case quantified along a line from the intracardial center (labeled "0") and extending to the ventricular lining (labeled "1") for various intertrabecular regions and trabeculae heights. The velocity magnitude strictly decreases from the center until around the neighboring trabeculae heights. The velocity magnitude then increases towards the center of the intratrabecular region, in some cases an order of magnitude, before dropping towards zero at the ventricle lining.
Interestingly, similar results are seen in the Re = 10 case, see Figure A6 in Appendix E. However, in the case of Re = 100, slightly different quantitative behavior is observed, see Figure 18. Note that the local minimum in velocity magnitude is still observed at the neighboring trabeculae height away from the ventricular wall, followed by an increase and then decrease as one moves closer to the wall. In contrast, as you measure away from the intracardial center, the velocity magnitude decreases (as before), but it also increases and then decreases before reaching the trabeculae height distance from the ventricle wall. This is due to the presence of a large intracardial vortex that forms, as illustrated in Figure 16. Each subfigure (A-E), corresponds to a different line between the intracardial center and ventricular lining in which the magnitude of velocity is measured. The velocity magnitude decreases from the center and then increases before decreasing again as one approaches a distance from the wall that is equal to the neighboring trabeculae heights. As one moves between the trabeculae, the velocity magnitude again increases, in some cases an order of magnitude, before dropping towards zero at the ventricle lining.

Pulsatile Flow through Idealized Trabeculated Chambers
Next, we consider the same idealized trabeculated ventricle but use a pulsatile inflow condition, as described in Appendix A. The pulsation frequency is given by a dimensionless frequency close to that reported for a 4 dpf embryonic zebrafish (f = 1.0). The Re was set to 0.  Next, we considered the horizontal velocity extending from the intracardial center to the intertrabecular region directly above it in the case of Re = 1. These results are shown in Figure 21. Flow velocities are at a minimum near the ventricular lining, with the magnitude of velocities decreasing at a slightly accelerated rate for larger trabeculae. On the other hand, at the intracardial center, the horizontal velocity increases with larger trabeculae in the middle of a pulsation cycle. As the pulsation cycle resides, it is clear that the flow direction also changes, giving rise to an intracardial vortex, as detailed in Figure 20.  Figure 22 shows the magnitude of velocity as as one moves along a line drawn from the top of the ventricle to the base of the middle intratrabecular region for five trabeculae heights five times during the pulse cycle. Note that, during a pulsation cycle, the snapshots taken during the middle of the pulse are similar to that of the steady inflow case for Re = 1 shown in Figure 16. For larger trabeculae heights, the horizontal velocity reaches a minimum at distance that is about a trabecular height away from the ventricle wall. As one moves between the trabeculae, the magnitude of the velocity increases and then approaches as one nears the wall. These velocity profiles confirm the formation of vortices as shown in Figure 20.    The spin direction is determined from the laminar flow, given in blue, moving left to right. When vortices first form in the intertrabecular regions, the vortices form moving counterclockwise. However, when a large intracardial vortex forms before intertrabecular vorticies, its direction is counterclockwise, which forces any intertrabecular vortices that form to move clockwise.
We also considered the horizontal velocity extending from the intracardial center to the intertrabecular region directly above for Re = 10, as shown in Figure 24. Similarly quantitative behavior is seen as in the Re = 1 case, see Figure 21; there is still significantly less flow in the intertrabecular region. However, the flow velocities are significantly different in the cavity below the intertrabecular region because of the presence of an intracardial vortex. It is clear that there is flow reversal given the differences in sign of the horizontal velocity. The results of the inertial dominated case, Re = 100, are shown in Figure 25. In all cases, a large intracardial vortex that fills the entire chamber is observed at the end of the pulse and the beginning of the next pulse. As the flow accelerates, the intracardial vortex is pushed downstream, and another intracardial vortex begins to form (t = 0.4T − 0.5T). One or more oppositely spinning vortices form between the trabeculae or between the two counterclockwise spinning intracardial vortices when t = 0.5T. The upstream intracardial vortex combines with the original intracardial vortex such that one large intracardial vortex is observed around t = 0.7T. When this occurs, the oppositely spinning vortices are annihilated. For h T b V ≥ 0.08, oppositely spinning intertrabecular vortices reappear at the end of the pulse.
Due to the formation of large intracardial vortices, the horizontal velocity changes sign when measured from the center of the cavity and proceeds directly upward toward the ventricle lining (see Figure 26). Similar quantitative behavior is seen as in the Re = 10 case. A substantial difference is the presence of a intracardial vortex that remains largely throughout the pulsation cycle. Moreover, similar to the other cases of Re = 1, 10, the velocity is significantly decreased within the intertrabecular region, even for Re = 100.

Conclusions
Two-dimensional immersed boundary simulations were used to solve for the flow fields within both biologically realistic geometries (3 dpf, 80 hpf, two 5 dpf WT zebrafish and a 7 dpf ErbB2-inhibited zebrafish from [21]) and idealized models of trabeculated ventricles. Specifically, we investigated the intracardial and intertrabecular fluid dynamics searching for possible vortex formation and spatially-varying velocity gradients.
The primary result of our parametric study using simplified models of the trabeculated embryonic heart is that small changes in morphology and the effective viscosity of the blood can result in significant changes in bulk flow patterns as well as the magnitude and direction of shear felt by the endocardium. This presents an interesting challenge since each of these parameters is continuously changing during development. The results also highlight the importance of obtaining geometric descriptions of the developing heart with high spatial and temporal resolution. Small differences in morphology between individuals may also result in fundamentally different flow signals, underscoring the need to quantify variation. Furthermore, the effective viscosity of the blood is indirectly proportional to Re, and small changes in Re may also result in drastic changes in flow direction and shear. The effective viscosity of the blood during development likely changes with hematocrit, which also changes from about 10-50% during development, and measurements of embryonic blood rheology are sparse [5].
This work focused specifically on the presence or absence of vortices given their significance to both the magnitude and direction of flow as well as the mixing patterns within the ventricle. When an intracardial vortex forms, the direction of the flow changes. When an intracardial vortex forms in unsteady flow, the direction of flow can change during the beat cycle, and the stagnation point moves along the cardiac wall. Since endothelial cells are known to sense and respond to changes in both magnitude and direction of flow, the formation and motion of these vortices could be important epigenetic signals.
The simulations revealed unexpected complexities in vortex dynamics as bulk flow moves through the chamber. In the cases of biologically realistic geometries, no large intracardial vortices developed at the biologically realistic fluid scale (Re ∼ 1) in the steady inflow case; in the pulsatile case, a large intracardial vortex formed when inflow velocities were minimal between pulses. However, for both steady and pulsatile inflow, intertrabecular vortices formed in pronounced trabeculated regions ( Figures 5 and 7). While similar behavior was observed in both the steady and pulsatile inflow cases, more vortices formed in the unsteady flow cases in regions with less pronounced trabeculation. These results were consistent with the idealized trabeculation model as well (Figures 16 and 20). When vortices form in intertrabecular regions, the flow changes direction in those areas compared to the direction of bulk flow in the chamber. In cases where not all intertrabecular spaces have a vortex, the flow between different trabeculae will move in different directions.
In addition, as expected, the velocity tapers off when measuring flow speeds close to the ventricular lining. In most cases, the flow tapers off approximately three-orders of magnitude before reaching an intertrabecular region. Interestingly, in these regions where there is pronounced vortex formation, the velocity increases and then decreases as it nears the wall (for the pulsatile cases with realistic geometries, see  Note that the corresponding steady inflow cases show similar trends). In the 5 dpf WT (I) case, the velocity may increase an order of magnitude during portions of a pulsation cycle (see Figure 10), while, in other cases, the increase is modest, such as in the 3 dpf WT case (see Figure 11). These results are quantitatively similar to those of the idealized model case at Re = 1 for both pulsatile and steady inflow. The presence of trabeculae appears to control the fluid velocities in this regions, which helps govern the amount of shear-stress felt at the endothelial layer.
It is almost certainly the case that the exact flow patterns and magnitudes will not be exactly the same for realistic geometries with moving walls. Some similar features are observed between our simulations and previous work using idealized geometries with moving walls [40] as well as realistic geometries with moving walls [19,20]. In all cases, there is a severe reduction in flow magnitude in the intratrabecular regions and the presence of intratrabecular vortices exhibiting changes in flow direction and magnitude as a function of intratrabecular depth. The consequences of the reduction in flow magnitude and the presence of vortices are that lower shear stresses are felt at the endothelial wall than would be the case without trabeculae and that the direction of shear changes both spatially and temporally. Direct comparison of our results to work with a moving boundary and idealized geometry [40] shows that our simplified models do not capture the exact spatial and temporal pattern of flow reversals and shear, but they do demonstrate that these features exist and are sensitive to perturbations in kinematics and morphology.
The idealized model cases expanded the study by allowing us to easily manipulate the system to understand its sensitivity to vortex formation, due to its complex geometry, fluid scale, and inflow characteristics. In this vein, we increased the span of our larger fluid scale from Re = 0.1 to Re = 100, rather than simply the biological case of Re = 1 for a 4 dpf embryonic zebrafish heart. We also investigated equally sized idealized trabeculae in a chamber, and varied the heights from no trabeculae to trabeculae twice as large as biologically relevant according to Figure 3a, and included both pulsatile and steady inflow simulations to parse effects of unsteady flows on vortex formation.
Moreover, diversity of trabeculated hearts across the animal kingdom are far and wide, and may cover a large spectrum of length and fluid (Re) scales. Even some invertebrate hearts contain ventricular trabeculation. In those invertebrates, their heart's morphology resembles that of lower vertebrates with sedentary lifestyles [50][51][52]. Some anthropods and mollusks contain ventiruclar trabeculae, such as blue crabs (Callinectes sapidus) [45], bar clams (Spisula solidissima) [53], oysters [54], snails [55], as well as octopus and squids [46]. However, any quantitative measurements detailing trabeculae morphology and flow measurements at varying time-points are unknown. It is possible that, in these hearts, the idealized simulations at higher Re or for larger trabecular heights are relevant to these invertebrate hearts.
The idealized models showed that a large intracardial vortex forms around Re ≈ 20 when steady flow is pushed through the chamber, while a similar sized vortex forms for Re = 10 when the flow is pulsatile. In general, pulsatile flow lowers the Re and trabeculae height needed to generate vortices. For both steady and unsteady flows as the trabeculae grow into the chamber, another bifurcation occurs in which small vortices form between each trabecula. Depending upon the Re and the morphology, the intertrabecular vortices can form without the presence of a large intracardial vortex (see Figure 16 for steady cases or Figures 19 or 20 for unsteady flow cases of Re = 0.1 or 1.0, respectively). In other cases, typically at higher Re, both the intracardial and intertrabecular vortices form, see Figures 23 or 25, for unsteady inflow for Re = 10 or 100, respectively. In all corresponding cases, the presence of large intracardial vortices changes the direction of the intertrabecular vortices. Note that, in the biologically relevant case of Re = 1 intracardial vortices do not form; this is consistent with the biologically accurate geometries as well. We do note that the exact patterns of intracardial and intratrabecular vortex formation may not hold for realistic geometries and pumping kinematics.
It is evident that there is a strongly coupled relationship between intracardial hemodynamics, genetic regulatory networks, and cardiac conduction. Besides contractions of the myocardial cells, which in turn drive blood flow, hemodynamics are directly involved in proper pacemaker and cardiac conduction tissue formation [56]. In addition, shear stresses are found to govern the conduction velocity distribution of action potentials within the myocardium [22]. Any changes in the emrbyonic heart's conduction properties will also affect the intracardial shear stresses, pressures, patterns of cyclic strains, and advection of morphogens. It is indeed a chicken and the egg scenario, especially when considering the first experiments that saw the importance of fluid dynamics in heart morphogenesis were performed in chicken embryos [57]. Dedicated initiatives to decipher exact cellular signalling pathways and genetic regulatory networks may be able to help further parse the causes of cardiac dysfunction. While CFD provides a robust framework to extract cardiac flow information (fluid velocities, shear stress distributions, pressures, etc.), coupling this data into a multi-scale cellular model is imperative to better understand the causes of many congenital heart diseases.
Unfortunately, the exact mechanisms of mechanotransduction are not yet clearly understood [58,59]. Biochemical signals are thought to propagate throughout a pipeline of epigenetic signaling mechanisms, which may regulate of gene expression, cellular differentiation, proliferation, and migration [60]. In vitro studies have discovered that endothelial cells detect shear stresses as low as 0.2 dyn/cm 2 , resulting in up or down regulation of gene expressions [61]. Shear stresses around ∼8-15 dyn/cm 2 are known to cause cytoskeletal rearrangement [62]. The aforementioned shear stresses reported in this paper and other CFD studies [20] are well in the range of those measured within embryonic hearts, ∼2 dyn/cm 2 and ∼75 dyn/cm 2 at approximately 1.5 and 4.5 dpf, respectively [3]. Mapping out the connection between fluid dynamics, the resulting ventricular stresses, electrophysiology, and the mechanical regulation of developmental regulatory networks are paramount to moving towards a more holistic understanding of heart development.
The simplified models used in this study permitted the exploration of a wide parameter space that includes the diversity of trabeculated hearts in vertebrate embryos and also in some invertebrates. There are obviously limitations in using such simplified models. Flow through the complex trabeculae of the ventricle is inherently three-dimensional. The flows generated by a beating ventricle will be different from pulsatile flows driven through a fixed geometry. Accordingly, our results should motivate the further development of more sophisticated three-dimensional models with complex geometry and moving boundaries. In the development of such models, our results illustrate how critical it is to obtain high resolution geometries and spatially and temporally resolved kinematics. Our results also point towards the consequences of variation in both geometry and blood rheology, which suggest that flow patterns may change drastically through development, between individuals, and across species. range, defined here as 0.01 < Re < 1000 (see, for example, [64][65][66][67][68][69]). However, IB is capable of solving fully coupled fluid-structure interaction systems, here we only use it to solve fluid flow through complex model geometry. The model consists of stiff boundaries that are immersed within an incompressible fluid of dynamic viscosity, µ, and density, ρ. The fluid motion is described using the full 2D Navier-Stokes equations in Eulerian form, given a where u(x, t) = (u(x, t), v(x, t)) is the fluid velocity, p(x, t) is the pressure, F(x, t) is the force per unit volume (area in 2D) applied to the fluid by the immersed boundary, i.e., the ventricle geometry. The independent variables are the position, x = (x, y), and time, t. Equation (A1) is equivalent to the conservation of momentum for a fluid, while Equation (A2) is a condition mandating that the fluid is incompressible. The interaction equations between the fluid and the immersed structure are given by where X(r, t) gives the Cartesian coordinates at time t of the material point labeled by Lagrangian parameter r, f(r, t) is the force per unit area imposed onto the fluid by elastic deformations in the boundary, as a function of the Lagrangian position, r, and time, t. Equation (A3) applies a force from the immersed boundary to the fluid grid through a delta-kernel integral transformation. Equation (A4) sets the velocity of the boundary equal to the local fluid velocity. The force equations are specific to the application. In a simple case where a preferred motion or position is enforced, boundary points are tethered to target points via springs. The equation describing the force applied to the fluid by the boundary in Lagrangian coordinates is given by f(r, t) and is explicitly written as f(r, t) = k target (Y(r, t) − X(r, t)) , where k target is the stiffness coefficient, and Y(r, t) is the prescribed Lagrangian position of the target structure. In all simulations, the immersed structure was held nearly rigid by applying a force proportional to the distance between the location of the actual boundary and the preferred position. The deviation between the actual and preferred positions can be controlled with the variable k target . The fluid flow is driven through the immersed boundary using either pulsatile parabolic inflows or a linear ramp to steady parabolic inflow at the location of the AV canal. The equations describing the specific inflow boundary conditions are given in Table A1. A partial Neumann outflow condition is enforced in the direction of flow at the outlet. This outflow condition is given as where u and v are the xand y-components of the fluid velocity, respectively, and ∂u ∂n is the directional derivative of the x-component of the velocity taken in the direction normal to the boundary of the fluid domain. Table A1. Inflow boundary conditions for both simulations, one pertaining to parabolic steady inflow and the other corresponding to a parabolic pulsatile inflow. The parameters used for the boundary conditions are f , the non-dimensional frequency, which is matched to the zebrafish heart at 96 hpf, and V in , the maximum inflow velocity.

Case
Inflow BC

Steady Inflow
We used an adaptive and parallelized version of the immersed boundary method, IBAMR [70,71], to perform the simulations involving the highly idealized trabecular models, for both the steady inflow and pulsatile inflow cases. IBAMR is a C++ framework that provides discretization and solver infrastructure for partial differential equations on block-structured locally refined Eulerian grids [72,73] and on Lagrangian (structural) meshes. IBAMR also includes infrastructure for coupling Eulerian and Lagrangian representations.
The Eulerian grid on which the Navier-Stokes equations were solved was locally refined near the immersed boundaries and regions of vorticity with a threshold of |ω| > 0.05. This Cartesian grid was organized as a hierarchy of four nested grid levels, and the finest grid was assigned a spatial step size of dx = D/1024, where D is the length of the domain. The ratio of the spatial step size on each grid relative to the next coarsest grid was 1:4. The temporal resolution was varied to ensure stability. Each Lagrangian point of the immersed structure was chosen to be D 2048 apart (twice the resolution of the finest fluid grid).

Appendix B. Methods: Biologically Realistic Geometries
We present two more computational geometries we will analyze for WT zebrafish embryos at 80 hpf and 5 dpf (II). Note that the 80 hpf is close to the 3 dpf time-point and the 5 dpf (II) zebrafish illustrates natural variation in WT zebrafish from the embryo given in Figure 2. Figure A1. Extracting biologically realistic geometries of trabeculated geometries at different time-points from [21]. Geometries are given for a WT zebrafish at 80 hpf and 5 dpf (II).

Appendix C. Results: Pulsatile Inflow for Biologically Realistic Geometries Results
Here, we give more velocity measurements extending from the AV canal center to intertrabecular regions along the ventricular canal for cases of a 3 and 5 dpf WT (I) and 7 dpf ErbB2 inhibited zebrafish, Figures A2-A4, respectively. These figures are given to compare to their cases of steady inflow. Note that, in all cases, the velocity profiles taken at 50% of the pulsation cycle looks qualitatively identical to its corresponding steady inflow case. Note that this data was presented in Figures 10-12 as semi-log plots to illustrate the extent of the velocity decay as measured near the ventricular lining.

Appendix D. Results: Computing WSS, TAWSS, and OSI
WSS was computed as in [74], e.g., letting σ = −pI + µ ∇u + (∇u) T denote the fluid stress, the WSS, τ can be given by wheret andn is the anti-clockwise tangential and normal vector, respectively, [σ] denotes the jump in the stress across the immersed boundary, and |∂X/∂s| gives the interfacial stretching. TAWSS is calculated by temporally averaging the WSS, τ, e.g., where T is the period of a cardiac cycle, as before. OSI was computed as in [19]; e.g., We accounted for noise in the force densities along the immersed boundary in two ways. First, the WSS and TAWSS were spatially averaged either across an entire trabeculae or intertrabecular region and the entire ventricle; we do not present shear-stress measurements at individual Lagrangian points but only at subsets of the Lagrangian domain. Second, the resulting OSI was spatially-averaged on each Lagrangian node using the closest five neighboring Lagrangian nodes.
In this appendix, we present magnitude of velocity results in the idealized geometry case for Re = 10 (see Figure A6 below). Figure A6 shows similar qualitative trends to that of the Re = 1 case in Figure 17 of Section 3.3. However, differences arise in the actual quantitative magnitude of velocity measurements, where, in the Re = 10 case, the velocity magnitudes are generally slightly greater. Figure A6. Measurements of the non-dimensional magnitude of velocity are provided for the Re = 10 case quantified along a line from the intracardial center (labeled "0") and extending to the ventricular lining (labeled "1") for various intertrabecular regions and trabeculae heights. The velocity strictly decreases from the center until around the neighboring trabeculae heights, in which the velocity increases, in some cases, an order of magnitude, before dropping towards zero at the ventricle lining.