1 Introduction

The morphogenesis and growth of the developing heart are influenced by hemodynamic forces, which are exerted on the heart walls by the flow of blood [1]. Experiments on embryos from chicks, zebra fish, and mice during early development have shown that deviations from normal blood flow in the embryo heart lead to malformations [2–5] that resemble heart defects in humans with congenital heart disease (CHD). CHD affects about 1% of all human live babies in the United States [6, 7].

Both in vitro and in vivo studies have shown that the cells that form the walls of the heart (and blood vessels) respond to alterations in hemodynamic forces by generating a cascade of signaling and gene expression events that ultimately lead to changes in heart morphology and growth [8–14]. Hemodynamic forces can be divided into: (1) blood pressure, which acts perpendicularly to the wall, and (2) wall shear stress, which acts tangentially to the wall. Wall shear stress originates from viscous (friction) forces due to blood movement near the wall and is proportional to the gradient of blood velocities in the direction normal to the wall. In order to predict the effect of abnormal blood flow on heart development, a better understanding of the mechanisms by which hemodynamic forces affect heart development is needed. A step towards achieving this goal is to characterize the dynamic distribution of hemodynamic forces (blood pressure and wall shear stress) acting on the walls of the developing heart over the cardiac cycle.

Inside the developing chick heart, blood pressure can be measured in vivo [e.g., 15, 16], and wall shear stress can be calculated from blood velocity profiles measured in vivo near the heart wall. Blood velocities in the hearts of chick embryos [17] and zebra fish [4] have been measured using micro-particle imaging techniques. However, such measurements present several challenges due to a combination of the small length-scales of the developing heart (<2 mm) and the continuous movement of the heart walls. Such difficulties, and their effect in the accuracy of the measurements, affect the calculation of wall shear stress in embryonic hearts [18]. Finite element models (FEM) have also been used to predict wall shear stress on the chick embryonic heart [14, 19, 20]. However, these previous FEMs were static, neglecting the dynamic effects of the heart wall motion and pulsatile blood pressure on wall shear stress.

Our objective is to determine how changes in the heart geometry and wall motion affect wall shear stress during the cardiac cycle, focusing on models of the outflow tract (OFT) of the chick heart during an early developmental stage (Hamburger Hamilton stage 18 [21]; HH18). To this end, we used dynamic FEMs, which differ in the geometries of the OFT, and that incorporated the effects of dynamic wall motion and pulsatile blood pressure on blood flow. Geometric parameters and motion patterns of the chick OFT walls used in the FEMs were estimated from high-resolution optical coherence tomography (OCT) images. We chose the chick for our animal model because (1) chick embryos are easy to access for imaging and measuring; and (2) in early stages of development, the chick heart resembles the human heart [22]. We chose to focus on the OFT of the chick heart at HH18 because at this stage the OFT has a relatively simple geometry suitable for biomechanical modeling [23, 24] and because the OFT is very sensitive to changes in hemodynamic conditions [3]. Quantification of the changes in wall shear stress over a cardiac cycle in the OFT provides a step towards elucidating the role of hemodynamic forces on heart development.

2 Biological problem

During development, the morphology of the chick embryonic heart changes from a tubular structure into a four-chambered heart [23]. The chick heart starts beating at HH10 (~36 h of incubation) [23, 25]. At HH18 (~68 h incubation), the heart of the chick embryo consists of a looped tube that pumps blood presumably via a peristaltic-like contraction motion [16, 26]. Heart septation and chamber formation start after HH21 (~84 h of incubation).

At HH18, the chick heart consists of contiguous segments: the sinus venosus (inflow tract), the primitive atrium, the atrioventricular (AV) canal, the primitive ventricle, and the OFT [23, 27]. Although the heart has no valves at HH18, the AV canal and the OFT have cardiac cushions, which are protrusions of the heart walls toward the lumen (domain where the blood flows) that are located where heart valves will form [23, 24, 28]. These cardiac cushions presumably increase the pumping efficiency of the developing heart. At HH18, the OFT is a slightly curved tube with an average external diameter of ~430 μm [29] and an approximate length of ~600 μm [27]. Typical heart rates of chick embryos at HH18 are 2.2–2.4 beats per s [16, 30]; thus, the period of the cardiac cycle, T is ~0.45 s. During ventricular systole (about 1/3 of the cardiac cycle [16]), when the ventricle is contracting and ejecting blood into the chick arterial system, the OFT is fully open, allowing the flow of blood from the ventricle to the aortic sac (Fig. 1a, c). However, during ventricular diastole (about 2/3 of the cardiac cycle), when the ventricle is filling with blood mainly from the atrium, the OFT walls contract (Fig. 1b, d), limiting reverse blood flow (backflow).

Fig. 1
figure 1

Sections of HH18 chick heart OFT obtained from OCT images. The figure shows OCT images of the OFT during the cardiac cycle: a, c during ventricular systole, when the OFT is most expanded, a is a longitudinal section and c is a cross-section; and b, d during ventricular diastole, when the OFT is most constricted, b is a longitudinal section and d is a cross-section. Point P marks the approximate location where velocities were acquired with Doppler OCT (see Fig 5.b). The white-dotted line corresponds to the direction of the incident light beam of Doppler OCT; the arrows indicate the direction of blood flow. Scale bar = 100 μm. L Lumen, M Myocardium; CJ Cardiac jelly

Full size image

The walls of the OFT are composed of three concentric layers [24]: endothelium, cardiac jelly, and myocardium. The endothelium, a single layer of endothelial cells (ECs), lines the internal part of the OFT wall and thus, it is in direct contact with blood flow. The cardiac jelly is comprised of an amorphous extracellular matrix that constitutes the bulk of the wall, including the cardiac cushions. The myocardium, consisting of a layer of primitive myocardial cells (MCs), is in the outer part of the OFT and actively contracts to limit backflow during ventricular diastole (note that in contrast to the OFT myocardium, the myocardium in the ventricle contracts to eject blood during ventricular systole). Previous research suggested that MCs respond mainly to changes in blood pressure [11–13], whereas ECs respond mainly to changes in wall shear stress [8–10, 14]. In vitro studies of ECs further suggested that the response of ECs depends on both changes in wall shear stress that occur over time (temporal variation) and the spatial distributions of wall shear stress [8–10, 31, 32]. To better understand the response of ECs to blood flow in vivo, we need a more comprehensive characterization of wall shear stress acting on the cardiac walls of the chick heart.

3 Methods

To determine the influence of wall geometry and motion on wall shear stress in the developing chick heart, we used dynamic FEMs of the chick OFT at HH18. Our OFT models were based on in vivo 2D images of the chick heart that captured the motion of the OFT during the cardiac cycle. Simulations of these models gave the temporal variations and spatial distributions of wall shear stress on the OFT during the cardiac cycle.

Previously, we developed a dynamic, 3D image-based FEM of the OFT based on the assumption that the OFT has circular lumen cross-sections [33]. However, at HH18 the OFT has cardiac cushions that render the lumen cross-section non-circular [28] (see Fig. 1c, d). To study the effect of cardiac cushions on the distribution of wall shear stress, we developed FEMs of the OFT with and without cardiac cushions.

3.1 Heart morphology and blood flow imaging

An OCT imaging system with a spatial resolution of ~10 μm was used to acquire 2D morphological images of the OFT of chick embryos at HH18 (n = 4). Details of the OCT system, which was based on a spectral domain configuration, were reported previously [34–36].

To image the OFT of the chick heart, fertilized white leghorn eggs were incubated with blunt end up at 102°F and 85 to 87% humidity, in a horizontal rotation incubator (No. 1536E GQF Mfg. Co., Savannah, GA) for 3 days. To access the embryonic heart, the egg shell was opened, and the membrane that overlays the chick heart was removed. The egg was then placed on a custom-made stage under the OCT probe and the embryo was gently positioned so that the OFT could be easily imaged. For each embryo, OCT acquired 2D images of longitudinal and transverse cross-sections of the OFT (Fig. 1) over a period of ~5 s, at a rate of 20 images per second. Thus, we captured ~10 frames per cardiac cycle. The images revealed the contraction and expansion of the OFT wall during the cardiac cycle and the presence of the cardiac cushions, which greatly reduced the cross-sectional area of the lumen during contraction.

In addition, velocity of blood flow inside the OFT was measured by the OCT system configured in Doppler mode [36]. Velocity data were acquired at a point located at the approximate center of the OFT (point P in Fig. 1a) with a time resolution of 0.1 ms. Measured velocity data corresponded to the projection of the 3D blood velocity vector in the direction of the incident OCT light beam (dotted line in Fig. 1a). Although measured velocity data did not provide an accurate description of the 3D velocity field inside the OFT, they revealed temporal variations of blood velocity over the cardiac cycle.

3.2 Mathematical model and finite element discretization

3.2.1 Mathematical model

Blood flow inside the chick OFT was modeled as an incompressible Newtonian fluid (i.e. using the Navier–Stokes equations). Because the OFT walls contract and expand during the cardiac cycle, the lumen geometry (fluid domain) changes continuously with time. Dynamic changes in the geometry of the OFT lumen affect blood velocities and therefore need to be incorporated into the equations of blood flow. This was done through an arbitrary Lagrangian–Eulerian (ALE) formulation of motion (e.g. [37, 38]), in which equations are expressed in terms of a moving reference frame (in the finite element implementation, the moving reference frame is the deforming finite element mesh). Using the ALE formulation, the Navier–Stokes equations (in Cartesian coordinates and using indicial notation) are

$$ \rho \left[ {\frac{{\delta v_{i} }}{\delta t} + (v_{j} - \hat{v}_{j} )v_{i,j} } \right] = \tau_{ij,j}^{F} $$
(1)
$$ v_{i,i} = 0 $$
(2)

where the indices i and j indicate components (in the x, y and z directions) and “,” indicates differentiation; ρ is density; v i indicates the (i)th component of the fluid velocity vector (with respect to a frame fixed in space); τ ij F is the (i,j)th component of the fluid stress tensor; δv i /δt is the time derivative of v i with respect to the moving reference frame (as measured by an observer moving with the frame); and \( \hat{v}_{i} \) is the (i)th component of the velocity of the moving reference frame (with respect to a frame fixed in space).

The constitutive relations for a Newtonian fluid are

$$ \tau_{ij}^{F} = - p\delta_{ij} + 2\mu e_{ij} $$
(3)
$$ e_{ij} = \frac{1}{2}\left( {v_{i,j} + v_{j,i} } \right) $$
(4)

where p is hydrostatic pressure, μ is viscosity, and δ ij is the Kronecker delta.

In this study, the walls of the OFT were only included in the models to simulate the deformations of the lumen cross-section with time. Wall motion was assumed to be quasi-static (transient terms were neglected), with governing equations of motion given by:

$$ \varvec{\tau}_{ij,j}^{S} = 0 $$
(5)

where \( \tau_{{_{ij} }}^{S} \) is the (i,j)th component of the wall Cauchy stress tensor. The OFT walls were assumed to be elastic. For an elastic material, the constitutive relations are:

$$ \tau_{ij}^{S} = \frac{E\nu }{(1 + \nu )(1 - 2\nu )}\varepsilon_{kk} \delta_{ij} + \frac{E}{1 + \nu }\varepsilon_{ij} $$
(6)
$$ \varepsilon_{ij} = \frac{1}{2}\left( {u_{i,j} + u_{j,i} } \right) $$
(7)

where E is the material Young’s modulus (elastic modulus), ν is Poisson’s ratio, and u i is the (i)th component of the wall displacement vector.

The coupling between the blood flow and wall equations was accomplished by satisfying two conditions (e.g. [38]):

  1. i)

    equilibrium of forces at the interface,

$$ \tau_{ij}^{F} n_{j}^{I} = \tau_{ij}^{S} n_{j}^{I} $$
(8)

and

  1. ii)

    compatibility

$$ \hat{u}_{i}^{I} = u_{i}^{I} $$
(9)

In Eqs. (8) and (9), the superscript I indicates interface, n j I is a unit vector normal to the wall that points towards the lumen, \( \hat{u}_{i}^{I} \) is the displacement of the lumen boundary, and u I i is the displacement of the wall at the interface with the lumen.

In addition, a no-slip condition was imposed at the interface between the lumen and the wall,

$$ v_{i}^{I} = \dot{u}_{i}^{I} $$
(10)

The no-slip condition ensured that blood particles in contact with the wall moved with the wall (at the same velocity).

The force per unit interfacial area exerted by the flow of blood on the walls of the OFT, the stress vector t i , was obtained from the fluid stress tensor (τ ij F) evaluated at the lumen-wall interface (e.g. [39]):

$$ t_{i} = \tau_{ij}^{F} n_{j}^{I} $$
(11)

where τ F ij is evaluated at the interface. t i can be decomposed into a normal vector and a vector tangential to the wall. The normal stress vector \( t_{{n_{i} }} , \) is the projection of t i into \( n_{i}^{I} \), and its magnitude is approximately equal to the hydrostatic pressure, p; hence, the normal vector is:

$$ t_{{n_{i} }} = (t_{j} n_{j}^{I} )n_{i}^{I} \approx - p\,n_{i}^{I} $$
(12)

The tangential stress vector or “wall shear stress vector”, \( \tau_{{w_{i} }} , \) is the projection of t i into the plane of the lumen-wall interface,

$$ \tau_{{w_{i} }} = t_{i} - t_{{n_{i} }} $$
(13)

Volume flow rate, Q, is the volume of blood flow that passes through a lumen cross-section per unit time. If \( n_{i}^{A} \) is a unit vector normal to the plane of the cross section (pointing towards the OFT outlet), and A is the area of the lumen cross-section, then

$$ Q = \int\limits_{A} {v_{i} } n_{i}^{A} {\text{d}}A $$
(14)

Q is positive when blood flows from the ventricle to the aortic sac and negative during backflow.

3.2.2 Geometry of OFT models

To investigate the influence of cardiac cushions on the wall shear stress vector \( \tau_{{w_{i} }} \) (eq. 13), we generated three models of the OFT: (1) a “cylindrical model”, in which the lumen was modeled as a straight circular cylinder (the walls were not explicitly modeled); (2) a “cushion model” (Fig. 2a, b) that included the OFT lumen and cardiac cushions, and (3) a “jelly model” (Fig. 2c, d) that included the OFT lumen, cardiac cushions, and an additional layer of wall (to simulate additional cardiac jelly material in contact with the myocardium).

Fig. 2
figure 2

Two FEMs of the OFT. Left: Cushion model: a reference cross-section and b FEM discretization. Right: jelly model: c reference cross-section and d FEM discretization. (Cylindrical model not shown.) In a and c, the inner part corresponds to the lumen and the outer part, to the wall. In b and d, the lumen-wall interface is marked with a thick line. Dimensions are: R = 0.1875 mm, L 1  = 0.077665 mm, L = 0.5 mm, h = 0.025 mm, and r = 0.03 mm. Points A and B are representative points where WSS was analyzed in detail

Full size image

In all three models, the external surface of the OFT was modeled as a straight circular cylinder. The cylindrical model, given its symmetry, was simulated using a 2D axisymmetric FEM of blood flow. The symmetry of the cushion and jelly models allowed us to represent only one quadrant of each model, which were simulated using a 3D fluid structure interaction procedure (see Fig. 2). The cylindrical, cushion, and jelly models differed in the geometry of the lumen cross-section, and therefore simulations of these three FEMs were used to determine differences in wall shear stress due to OFT lumen geometry.

In analyzing wall shear stress on the OFT using the cushion and jelly models, we focused on two points (see Fig. 2a, c): A and B. Because of the symmetry of these models, \( \tau_{{w_{i} }} \) at points A and B has always the same direction (relative to the “cell” position) but changes in magnitude during the cardiac cycle. For the cylindrical model, points A and B cannot be distinguished, and wall shear stress is uniform at the lumen-wall interface. We assumed that the “magnitude” of the wall shear stress vector \( \tau_{{w_{i} }} , \) WSS was positive during forward flow and negative during backflow.

3.2.3 Boundary conditions on the OFT models

3.2.3.1 Motion of the OFT wall

To simulate the passive distension and active contraction of the OFT myocardium during the cardiac cycle, a radial displacement (see Fig. 3) was prescribed on the OFT external surface, which had cylindrical symmetry for all three models. The temporal variation and amplitude of the prescribed radial displacement were estimated from the OFT cross-sectional images acquired with OCT over time (e.g., Fig. 1c, d) and simplified as shown in Fig. 3. The amplitude of the prescribed radial displacement, D1 was 70 μm (see also [29]). The reference geometry of the cushion and jelly models (when prescribed radial displacement was zero) is shown in Fig. 2. Maximum OFT expansion corresponds to a radial displacement of 0.5 D 1, and maximum OFT contraction corresponds to a radial displacement of −0.5 D 1. A different amplitude, D 2 = 90 μm, was also used to assess the effect of wall motion amplitude on blood flow and WSS.

Fig. 3
figure 3

Blood pressure and myocardium displacements prescribed as boundary conditions on the OFT models. Top temporal variations of ventricular pressure (P v) prescribed at the inlet surface and pressure prescribed at the outlet surface (P a). Bottom radial displacement prescribed on the external surface of the OFT models with amplitude D 1. Note the different scales for blood pressure (on left vertical axis) and radial displacements (on right vertical axis)

Full size image

To further study the influence of wall motion on blood flow dynamics, we simulated two types of wall motion: (1) Simultaneous motion, in which prescribed radial displacements were the same along the OFT longitudinal direction; that is, all points on the OFT external surface moved the same amount in the radial direction, at the same time. Simultaneous motion was applied to all three FEMs. (2) Peristaltic motion, in which prescribed radial displacements were modeled as a displacement wave traveling along the OFT longitudinal direction, at a velocity of 7 mm/s [26]. In other words, peristaltic motion involves time lags in the wall motions of contiguous OFT cross-sections. Peristaltic motion was applied to the cushion model. The simulated wall motions were used to determine the effect of wall motion on Q and WSS.

In our three OFT models, the walls were not allowed to expand or contract in the longitudinal direction. This was accomplished by restricting the longitudinal motion of the OFT ends.

3.2.3.2 OFT lumen inlet and outlet conditions

For the boundary conditions on the OFT lumen inlet, we used published ventricular blood pressure data [16]. We prescribed a simplified pulsatile pressure wave (see Fig. 3) on the lumen inlet surface of the OFT models, as a normal stress vector or normal traction.

Since there were no available data for blood pressure at the outlet of our OFT models (close to the aortic sac), we used blood pressure data in the chick dorsal aorta [15] to estimate blood pressure at the outlet. Because the location of the dorsal aorta is relatively far downstream from the outlet of our OFT models, using dorsal aorta pressure directly would have resulted in an over-estimation of the pressure difference between inlet and outlet, ΔP. To overcome this difficulty, we prescribed a simplified pulsatile pressure on the surface of the lumen outlet that preserved the temporal variation of blood pressure in the dorsal aorta, but in which maximum and minimum pressure values were increased to reduce ΔP. The increase in maximum and minimum pressure values (from dorsal aorta measurements) was determined by assuming ΔP to be equal to 0 at the start of ventricular ejection, and 33 Pa at ventricular systole, about the same mean ΔP measured in the AV canal during ventricular filling [16]. This later assumption was made because the AV canal and OFT apparently have similar roles in regulating blood flow through the developing heart [25, 30]. Outlet blood pressures were prescribed as normal tractions.

Blood pressures at the inlet and outlet surfaces of our OFT model were the same for all FEMs. Therefore the temporal variation of ΔP was the same for all OFT models.

3.2.4 FEM implementation

In our OFT models, embryonic chick blood was assumed to be a viscous, incompressible Newtonian fluid with a density of ρ = 1,060 kg/m3 and a viscosity of μ = 3 × 10−3 kg/m s [19, 33], and the flow of blood was assumed to be laminar [40]. The OFT walls in the cushion and jelly models were assumed to be an almost incompressible elastic material with a Poisson’s ratio of ν = 0.49 and a Young’s modulus of E = 1,000 Pa.

The lumen of our models was discretized using flow-condition-based interpolation (FCBI) elements [41], with the mesh near the wall slightly refined to better capture velocity variations near the wall and therefore, to calculate wall shear stress with greater accuracy. Using FCBI elements to discretize our models is equivalent to using control volume methods to calculate blood flow [42]. The specific FEM discretizations used for each of our OFT models are summarized in Table 1.

Table 1 Finite element discretization of the OFT models
Full size table

Blood flow was assumed to be initially at rest and three cardiac cycles were simulated. Each cardiac cycle (T = 0.45 s) was discretized using 100 time steps, with each time step equal to 4.5 × 10−3 s. All of our OFT models were simulated using the FEM software ADINA (Watertown, MA) [41].

3.2.5 Convergence study

To ensure the accuracy of the results obtained, we performed a convergence study. We focused on the convergence of the blood-flow solution since, in this study; we were primarily interested in calculating WSS (wall shear stress). The walls of the chick OFT were simulated to account for the approximate deformation of the lumen cross-section over a cardiac cycle, but not to calculate strains and stresses in the wall accurately. For the blood flow convergence study, we compared volume flow rate Q and WSS at representative points (points A and B in Fig. 2a, c). Reported values correspond to those obtained at the mid-cross-section of the OFT models.

Convergence of the cylindrical model was assessed by comparing the FEM solution—using simultaneous wall motions and quasi-static conditions (in which transient terms in the fluid equations are neglected)—with the solution of the Hagen–Poiseuille flow under identical ΔP and model radius, R. Using the Hagen–Poiseuille solution [43], derived for a fully developed flow in a circular tube, volume flow rate (Q) is equal to:

$$ Q = \frac{\pi }{8\mu }\frac{\Updelta P}{L}R^{4} $$
(15)

and wall shear stress (WSS) is:

$$ {\text{WSS}} = \frac{\Updelta P}{2L}R $$
(16)

where L is the length of the cylinder in the longitudinal direction. The maximum differences between the Hagen-Poiseuille solution and the results of the cylindrical model were 0.05% for Q and 2.5% for WSS, suggesting that the FEM mesh that was used for the cylindrical model captured the behavior of the blood flow with accuracy.

To test the convergence of the cushion and jelly models, the results obtained with the fluid mesh specified in Table 1 and a mesh with ~eight times more fluid elements were compared for the cases in which the OFT had the maximum expansion and contraction (see Fig. 4). Maximal differences found between the meshes were 2.7% for Q and 3.2% for WSS, suggesting that the meshes shown in Table 1 and that were used in this study were sufficiently dense.

Fig. 4
figure 4

Comparison of OFT lumen cross-sectional areas, obtained from our three models when a OFT is most expanded, and b OFT is most constricted

Full size image

We also simulated the transient behavior of the cushion and jelly models using the mesh specified in Table 1 but decreasing the time-step to 2.25 × 10−3 s (from 4.5 × 10−3 s used in the simulations). Maximal differences between solutions obtained with the two time steps were 0.5% for Q and 3% for WSS. Therefore, the time step employed was sufficient to provide accurate solutions.

4 Results

We report results as a function of a non-dimensional time t/T, where t is time and T is the period of the cardiac cycle (T = 0.45 s in our simulations). For reference, in our models ventricular systole is between t/T = 0 and t/T = 0.38, and ventricular diastole from t/T = 0.38 to t/T = 1. We chose the cushion model as our basic OFT model, against which we compared the cylindrical and jelly models.

4.1 Importance of inertial effects on blood flow

Blood flow in the chick heart is characterized by low Reynolds (Re) and Womersley (Wo) numbers [17, 33]. Re—a measure of the ratio between inertial and viscous forces [43]—is defined as: Re = ρ V R/μ, where ρ is the fluid density, V is a characteristic velocity, R is radius (or a characteristic length), and μ is the fluid viscosity. A small Re therefore implies that flow inertial forces are negligible. Wo—a measure of viscous effects in oscillating flows [43]—is defined as: Wo = ρ Ω R 2/μ, where Ω is a characteristic frequency (in our case the frequency of the cardiac cycle—about 2 Hz). A small Wo indicates that velocity is in phase with the pressure gradient. Flow with small Re and Wo numbers is laminar (not turbulent).

For the cylindrical model, we compared the results of simulations under quasi-static conditions with those of a transient analysis. These comparisons showed that the temporal variation of the center velocity, Q and WSS, were very similar for both cases (results not shown). The maximal differences occurred during contraction of the OFT wall, when the transient solutions presented a small time lag of about 0.003 s (less than a time step). Differences between the solutions obtained with the two approaches when the OFT was most constricted (t/T = 0.5) were 4% for the center velocity, 3.5% for Q and 2.8% for WSS. These small differences were in agreement with the small Re (maximum ~27) and Wo (maximum ~0.24) obtained from the transient solution.

Small Re and Wo numbers were obtained in all our simulations of OFT models. Maximal Re = 17.9 and maximal Wo = 0.26—for the cushion and jelly models, the “radius” R in the definition of Re and Wo numbers was assumed to be the distance between the center of the model cross-section and point B, see Fig. 2a and c. Obtained Re and Wo from our OFT models were consistent with previous reports [17, 33], and confirmed that inside the OFT inertial effects are negligible.

4.2 Blood velocity and volume flow rate

Figure 5a shows the axial velocity of blood flow at the center of the OFT cross-section during the cardiac cycle calculated from the cushion model (with simultaneous and peristaltic wall motions), and Fig. 5b shows the velocity measured using Doppler OCT at point P in Fig. 1a. Since the measured velocity is the projection of the blood velocity vector in the direction of the incident light beam (the dotted line in Fig. 1a), we could compare only the shape of the velocity profiles over time but not the absolute values of measured and calculated velocities. The calculated velocity (especially for the cushion model with peristaltic wall motion) resembled the temporal features of the measured velocity. Both calculated and measured velocities showed a peak during ventricular systole and significant backflow during ventricular diastole. However, the ratio between the absolute values of the peak positive and negative velocities was smaller for the calculated velocity (~1.4 for simultaneous and ~2.0 for peristaltic wall motion) than that of the measured velocity (~3.5), suggesting that backflow in our models was overestimated. In our OFT models, blood velocity just before ventricular systole was negative, whereas the measured Doppler OCT data showed some oscillations above zero. Discrepancies between calculated and measured velocities can be attributed to the simplifying assumptions of the cushion OFT model and to experimental uncertainties during velocity measurement, e.g. point P was fixed in space and therefore its relative position inside the OFT changed due to rigid body motions of the beating chick heart; the same is true for the angle between blood flow and the direction of the incident light beam. Given that our objective was to determine the influence of wall geometry and motion on WSS, the discrepancies between calculated and measured velocities were not crucial for our analysis.

Fig. 5
figure 5

Blood velocities calculated and measured inside the OFT. a Calculated longitudinal velocities at the center of the mid-cross-section obtained from the cushion model with simultaneous and peristaltic wall motions, with prescribed wall motion amplitude D 1 (70 μm). b Blood velocities measured with Doppler OCT at a fixed point located at the approximate center of the OFT lumen (see point P in Fig. 1a) in the direction of incident OCT light beam (dotted line in Fig .1a)

Full size image

Figure 6 shows Q obtained from simulations of the cylindrical, cushion and jelly models. The temporal variation of Q obtained from the cushion model was very similar to that from the jelly model (see also Table 2). In contrast, the magnitude of Q obtained from simulations of the cylindrical model was about two times larger than Q from the cushion model during ventricular systole and about four times larger during peak backflow. Differences in Q among the models, might be due to differences in the lumen cross-sectional areas of the OFT models. Simulations of the cylindrical, cushion, and jelly models resulted in variations in lumen cross-sectional areas during the cardiac cycle (see Fig. 4). Differences in the lumen cross-sectional areas between the cushion and jelly models were relatively small (about 5% difference when the OFT was most expanded). The cylindrical model, because of the absence of cardiac cushions, had a significantly larger lumen cross-sectional area than the cushion and jelly models (about 33% difference from the cushion model when the OFT was most expanded).

Fig. 6
figure 6

Volume flow rates (Q) at the mid-cross-section of the OFT obtained from the OFT models. For cylindrical, cushion, and jelly models, a simultaneous wall displacement of amplitude D 1 was prescribed. For the cylindrical SA model (Sect. 4.2), prescribed wall displacements were such that matched temporal variations of the lumen cross-sectional area calculated from the cushion model

Full size image
Table 2 Volume flow rates (Q) and wall shear stress (WSS) obtained from simulations of the OFT models
Full size table

To determine whether the geometry of the OFT lumen also affects Q, we varied the radius of the cylindrical model over time so that its cross-sectional area matched the cross-sectional area of the cushion model (we called this resulting model the cylindrical “same area” (SA) model). Simulations of the cylindrical SA model showed that Q was greatly reduced, when compared to the original cylindrical model (see Fig. 6). However, absolute maximum and minimum values of Q obtained from the cylindrical SA model were still larger than those from the cushion and jelly models (see also Table 2). These results suggest that cardiac cushions minimize backflow primarily by reducing the OFT lumen cross-sectional area but also by changing the geometry of the lumen.

Changes in the motion of the OFT wall also resulted in different temporal variations of Q. Figure 7a shows differences in Q obtained from the cushion model when the prescribed wall motion was simultaneous and peristaltic. Peristaltic motion reduced the absolute value of Q during backflow (see Table 2) and smoothed the abrupt variations of Q (produced by the non-smooth ΔP prescribed). Increasing the amplitude of OFT prescribed displacement from D 1 (70 μm) to D 2 (90 μm), affected Q (Fig. 7b). During ventricular systole, Q was larger when the amplitude of the wall motion was D 2 than when it was D 1 (see Table 2) because of an increase in the area of the OFT lumen. During backflow, however, the differences in Q between simulations of the two displacement amplitudes were smaller than during forward-flow.

Fig. 7
figure 7

Volume flow rates (Q) obtained from the cushion model. a Comparison between simulations of simultaneous and peristaltic OFT wall motions, calculated using a prescribed radial displacement amplitude D 1 = 70 μm. b Comparison between simulations of simultaneous OFT wall motion, calculated with prescribed radial displacement amplitudes D 1 = 70 μm and D 2 = 90 μm

Full size image

4.3 Wall shear stress

Variations in WSS were influenced by the geometry of the OFT models. While our cylindrical models, due to their symmetry, presented a uniform WSS along the perimeter of the lumen cross-section (and even along the OFT if we assume a fully developed blood flow profile), our cushion and jelly models showed a non-uniform distribution of WSS, with maximal WSS at the cardiac cushion (Fig. 8).

Fig. 8
figure 8

Temporal variations of WSS in the cushion model at points A and B (Fig. 2a) in the mid-cross-section of OFT. Results correspond to the case of simultaneous wall motion with prescribed wall displacement amplitude D 1 = 70 μm

Full size image

To evaluate differences in WSS, for the cushion and jelly models two points were considered: (1) point A, located at the center of the cardiac cushion, at the lumen-wall interface; and (2) point B, located at the center of the lumen-wall interface where there is no cardiac cushion (Fig. 2a, c). The temporal variations of WSS for the cushion and jelly models at points A and B were very similar (see Fig. 9). Temporal variations of WSS in the cylindrical SA model differed from those in the cushion and jelly models at points A and B, see Fig. 9. These results showed that the geometry of the OFT affected the temporal variation and spatial distribution of WSS.

Fig. 9
figure 9

Comparison of WSS obtained using the OFT models. WSS at a point A and b point B. The cushion and jelly models were simulated with simultaneous wall motion, with prescribed wall displacement amplitude D 1 = 70 μm; the wall motion of the cylindrical SA was prescribed such that temporal variations of cross-sectional area matched those of the cushion model

Full size image

Changes in OFT wall motion also resulted in variations of WSS. Peristaltic wall motion in the cushion model smoothed the temporal variation of WSS at points A and B and introduced slight changes in the amplitude of oscillation of WSS compared to simultaneous wall motion applied to the cushion model (see Fig. 10a, c and Table 2). Increasing the amplitude of the prescribed radial displacement from D 1 to D 2, using simultaneous wall motion, changed the magnitude of WSS at points A and B (Fig. 10b, d; and Table 2) with a larger change observed at point A during ventricular diastole. Maximal WSS during systole and minimal WSS during diastole; however, did not change significantly, even though Q changed significantly with changes in wall motion amplitude (see Fig. 7b and Table 2). These results suggest that changes in geometry affect the temporal variation and spatial distribution of WSS more than do variations in wall motion.

Fig. 10
figure 10

WSS obtained from the cushion model. Left: WSS at point A. a simultaneous versus peristaltic motion, and b prescribed wall displacement amplitude D 1 (70 μm) versus D 2 (90 μm). Right: WSS at point B. c simultaneous versus peristaltic motion, and b prescribed wall displacement amplitude D 1 (70 μm) versus D 2 (90 μm)

Full size image

4.4 Oscillatory shear index

The oscillatory shear index (OSI), following Ku et al. [44], estimates the change in direction of the wall shear stress vector,\( \tau_{{w_{i} }} \) (defined by eq. 13) with respect to a “mean” direction. This change in direction occurs due to changes in the direction of blood flow. Due to the symmetry of our OFT models, blood flow at points A and B is either forward or backward (no cross-flow can occur). Thus, OSI was defined as:

$$ {\text{OSI}} = \frac{{\int_{0}^{T} {{\text{WSS}}^{*} {\text{d}}t} }}{{\int_{0}^{T} {\left| {\text{WSS}} \right|{\text{d}}t} }}, $$
(17)

with

$$ {\text{WSS}}^{*} = \left\{ {\begin{array}{*{20}c} {0\quad {\text{if}}\,{\text{sign}}({\text{WSS}}) = {\text{sign}}(\int_{0}^{T} {\text{WSS}} {\text{d}}t)} \\ {\left| {\text{WSS}} \right|\quad {\text{if}}\,{\text{sign(WSS}}) \ne {\text{sign}}(\int_{0}^{T} {\text{WSS}} {\text{d}}t)} \\ \end{array} } \right. $$

where T is the period of the cardiac cycle. Positive WSS corresponds to forward flow of blood and negative to backflow. OSI then is the fraction of the integral of |WSS| in which WSS has a sign (+ or −) opposite to the sign of the mean. Note that if WSS is constant, then OSI = 0, and when WSS oscillates in a sinusoidal fashion, OSI → 0.5.

In all cases (except for the cushion model with prescribed radial displacement amplitude D 2), mean WSS was negative (Table 2). Therefore, OSI gave the ratio of the integral in time of positive WSS (during forward flow) and the integral in time of |WSS| in a cardiac cycle. It also provided a means to estimate the ratio between the integral of |WSS| during forward and reverse flow. We found that (see Table 2): (1) OSI was smaller in the cylindrical model compared to the cushion and jelly models; (2) OSI was larger at point A (i.e. at the cushion) than at point B (cushion and jelly models); and (3) peristaltic motion had almost no influence on OSI.

5 Discussion

The cylindrical, cushion, and jelly models of the HH18 chick heart OFT presented here were developed to determine the influence of OFT geometry and wall motion on WSS. The OFT in the cylindrical model had a circular lumen cross-section, and the cushion and jelly models (Fig. 2) had non-circular OFT lumen cross-sections because of the presence of cardiac cushions. For each model, radial displacements were prescribed on the external surface of the models, to simulate myocardium contraction/distension and wall motion during the cardiac cycle. Also for each model, a drop in blood pressure, ΔP, that changed over time was prescribed between the inlet and outlet surfaces of the OFT lumen (by applying “pressure” boundary conditions to the inlet and outlet surfaces, Fig. 3) to simulate the driving force for blood flow within the OFT.

Our OFT models were based on three main assumptions: (1) the pressure prescribed at the lumen outlet was estimated: since measured pressure data at the aortic sac (the outlet of our models) was not available, the prescribed outlet pressure was estimated from measurements downstream of the OFT; (2) the geometry of the FEMs was simplified: the OFT was modeled as a straight tube (instead of a bended tube), with either circular cross-sections (cylindrical model) or with a pair of identical and symmetric cardiac cushions (cushion and jelly models); and (3) the motion of the OFT wall was simulated by prescribing a radial displacement: the external surface of our models remained circular at all times and the temporal variation of radial displacements was estimated (from ten sequential OCT images). Given these assumptions, our models of the OFT neglected several aspects of actual blood flow inside the OFT.

Regarding Assumption 1, the estimation of prescribed outlet pressure affected ΔP; and since ΔP provides an important driving force for blood flow through the OFT, calculated Q and WSS were also affected. Therefore, our models provided only estimates of Q and WSS. Furthermore, in our OFT models, pressure boundary conditions were prescribed and hence uncoupled from the OFT dynamics. In living chick embryos, because cardiovascular circulation forms a closed system, ΔP and the geometry of the OFT lumen and wall motions are coupled. Uncoupling these variables, however, allowed for less complex calculations and interpretation of results. Regarding Assumption 2, the FEMs neglected the influence of the OFT curvature on blood flow, which affects the local variation of WSS. In addition, actual cardiac cushions are not as symmetric as in our models (Fig. 1), and therefore our models neglected non-symmetric distributions of WSS. Regarding Assumption 3, the motion of the OFT wall was assumed to go from an open to a close state (Figs. 3, 4), which is likely an oversimplification of the actual wall motion. In addition, the external surface of our OFT models was assumed to remain circular at all times, and therefore the effects of the non-symmetric motion of the myocardium layer, observed in the OCT images, were neglected in our simulations. Given these assumptions, our OFT models represent a starting point to estimate variations in Q and WSS during the cardiac cycle and to determine how wall geometry and motion change WSS in the chick developing heart.

Our results showed that volume flow rate, Q, was generally in phase with the imposed ΔP, but was also influenced by the motion of the OFT wall. Small phase shifts between Q and ΔP occurred due to wall motion, especially when the OFT was quickly expanding or contracting. Given the small length scales of the developing heart (<2 mm), blood flow inside the heart is strongly affected by viscous forces. This is reflected by the small Re (<30) and Wo (<1) numbers obtained from the simulations. Thus, Q was expected to generally be in phase with the imposed ΔP, but with amplitude modulations and small time shifts due to the effect of OFT wall motion on blood flow. We found that the amplitude of the wall motion contributed to the amplitude of oscillation of Q (Fig. 7b), as expected, and that peristaltic motion (due to the time lags in the motion of contiguous cross-sections) altered the temporal variation of Q with respect to that of simultaneous wall motion (Fig. 7a). The simulations performed, indicated that peristaltic motion limits backflow more efficiently than a simultaneous motion. In our models this is because for the simulated peristaltic motion the inlet of the OFT started contracting 0.036 s (0.08 T) before contraction started in the simulated simultaneous wall motion case.

Since WSS is affected by Q inside the heart, WSS depends on ΔP and the motion of the OFT wall. Changes of WSS over time were influenced by the amplitude of the wall motion (D 1 or D 2) (see Fig. 10b, d). We found, however, that wall motion had a larger influence on Q than on WSS (compare Figs. 7b with 10b, d). For a cylindrical model, the Hagen–Poiseuille solution (eqs. 15 and 16) predicts that WSS is proportional to Q 1/4. Thus, while variations in Q affect WSS, changes in WSS are relatively smaller than changes in Q. This is consistent with results from our models (although the proportionality with Q 1/4 is no longer valid for the cushion and jelly models).

Our results showed that cardiac cushions limit backflow in the OFT by constricting the lumen area and affecting the distribution of blood flow velocities in the lumen of the OFT. When the OFT myocardium contracts and pressure inside the OFT decreases, the cardiac cushions bulge towards the OFT lumen (see Fig. 4), reducing the area of the OFT lumen cross-section. This decrease in area significantly reduces Q during backflow—a reduction that is evident when comparing the calculated Q in the cushion and jelly models with the calculated Q in the cylindrical model (Fig. 6). The presence of cardiac cushions also changes the geometry of the OFT lumen cross-section, which we found to play a role in reducing backflow. Simulations of the cylindrical model, in which the cross-sectional area changed over time like the cushion model (cylindrical SA) showed that the magnitude of Q in models with cardiac cushion was smaller than in the cylindrical SA model, especially during backflow (Fig. 6). These results indicate that Q is influenced not only by the cross-sectional area of the OFT lumen but also by the geometry of the lumen cross-section. Differences in Q between models with and without cardiac cushions were consistent with the increase in the area of the lumen-wall interface in models with cardiac cushions (relative to cylindrical SA models). A larger interface increases the viscous resistance to blood flow (on the walls of the OFT), and therefore for the same applied ΔP and cross-sectional area, models with larger interfacial area result in lower Q. Hence, the geometry of cardiac cushions is effective in reducing backflow by reducing cross-sectional area and increasing the resistance to blood flow.

The presence of cardiac cushions also affects the temporal variation and spatial distribution of WSS at the lumen-wall interface. Because cardiac cushions change the geometry of the lumen cross-section (from that of a cylindrical tube), they affect Q and the distribution of blood velocities within the lumen and hence WSS (see Figs. 8, 9). The smaller OSI found for the cylindrical SA model compared to that of the cardiac and jelly models, indicates that the cardiac cushions increase the ratio between the time integral of |WSS| during forward and backflow (mean WSS in our models was generally negative). We also found that for models with cardiac cushions, at point A (cardiac cushion) the ratio of the integral of |WSS| during forward and backflow is smaller than at point B (OSI at point A is larger than at point B). Whether endothelial cells (ECs) can differentiate a positive from a negative WSS, or whether ECs respond to the magnitude of WSS, the integral of WSS over time or some other related variable is not yet known. However, it is clear that the ECs that line the cardiac cushion (where valves will later form) are subjected to different mechanical stimuli than cells outside the cardiac cushions. The biological implications of these spatial variations in mechanical stimuli are not yet fully understood, but WSS likely plays an essential role on valve formation.

6 Conclusions

We developed OFT models of HH18 chick embryos to characterize patterns of WSS over the cardiac cycle and to determine whether these patterns are affected by OFT geometry and wall motion (simultaneous and peristaltic wall motions were considered and two amplitudes of wall motion were simulated). Our OFT models showed that even though backflow through the OFT is limited by the presence of cardiac cushions, WSS during backflow can have a larger absolute value than WSS during forward flow. WSS showed significant oscillation over time and was asymmetrically distributed, with larger absolute values at the cushions. Therefore, ECs located at the cardiac cushions are subjected to different WSS than are ECs located elsewhere. These differences in WSS may affect valve formation.