Computational Model for Early Cardiac Looping

Abstract

Looping is a vital event during early cardiac morphogenesis, as the initially straight heart tube bends and twists into a curved tube, laying out the basic pattern of the future four-chambered heart. Despite intensive study for almost a century, the biophysical mechanisms that drive this process are not well understood. To explore a recently proposed hypothesis for looping, we constructed a finite element model for the embryonic chick heart during the first phase of looping, called c-looping. The model includes the main structures of the early heart (heart tube, omphalomesenteric veins, and dorsal mesocardium), and the analysis features realistic three-dimensional geometry, nonlinear passive and active material properties, and anisotropic growth. As per our earlier hypothesis for c-looping, actin-based morpho-genetic processes (active cell shape change, cytoskeletal contraction, and cell migration) are simulated in specific regions of the model. The model correctly predicts the initial gross morphological shape changes of the heart, as well as distributions of morphogenetic stresses and strains measured in embryonic chick hearts. The model was tested further in studies that perturbed normal cardiac morphogenesis. The model, taken together with the new experimental data, supports our hypothesis for the mechanisms that drive early looping.

INTRODUCTION

Cardiac looping is a crucial morphogenetic event during early embryonic development. During looping, the relatively straight heart tube (HT) bends and twists into a curved tube normally directed toward the right side of the embryo (Fig. 1). Because looping represents the first major visible indication of left-right asymmetry in the vertebrate embryo, and since abnormal looping can lead to a number of congenital heart defects,6,24 the mechanisms that drive and regulate this process have received considerable attention. Recent work has provided new insight into genetic, molecular, and biophysical facets of looping. However, much remains unknown about this complex process. The present work deals with the biomechanical side of this problem.

FIGURE 1.

Embryonic chick heart morphology during c-looping (ventral view). (A) Embryo at stage 10⁻, the heart is enclosed in a box and approximate dimensions are shown in micrometers. Fluorescent labels placed on the ventral midline help visualize dextral rotation. (B) Stage 10⁺. (C) Stage 11; labels move to outer curvature. HT=heart tube, C=conotruncus, OV=omphalomesenteric vein.

During looping, the HT transforms first into a c-shaped tube (c-looping) and then into an s-shaped tube (s-looping), laying out the basic plan of the four-chambered heart.18 This paper focuses on c-looping, as the HT undergoes a combination of ventral bending, dextral (rightward) torsion, and a relatively small amount of dextral bending.18,30

Prior studies in our laboratory suggest that the bending and torsional components of c-looping, although coupled to some extent, are driven by two different sets of forces.28 In particular, we have postulated that ventral bending is driven mainly by active myocardial cell-shape changes, possibly due to actin polymerization.13 Torsion, on the other hand, is caused by external forces supplied caudally by the omphalomesenteric veins (OVs) and ventrally by the splanchnopleure (SPL).30,31 The feasibility of these mechanisms has been demonstrated qualitatively through experiments and relatively simple computer models.

The purpose of the present study was to examine these proposed looping mechanisms quantitatively and to determine how active morphogenetic forces (e.g., actin polymerization, cytoskeletal contraction, and growth) are integrated to produce the form of the heart during the initial stages of c-looping. To facilitate this work, a computational model was developed for the embryonic heart. The model is based on realistic three-dimensional geometry, material properties, and boundary conditions. Morphogenetic force distributions are specified according to published experimental data. In addition, new regional stress and strain measurements in embryonic chick hearts are presented and used with the model to determine actin polymerization, contraction, and growth parameters, to evaluate the accuracy of the model, and to test the validity of our hypotheses for looping. The result is a model that helps us understand the complex interplay of different mechanical forces in the looping problem.

BIOLOGICAL BACKGROUND

The heart is the first functioning organ in the embryo. Because development of the chick heart essentially parallels that of the human heart, the chick embryo is a popular experimental model for studying cardiac morphogenesis. The system of Hamburger and Hamilton5 divides the 21-day incubation period of the chick embryo into 46 stages. Heart development begins at about stage 9 (29–33 h), as paired primordia fuse to form a straight, almost symmetric tube.

Looping starts at about stage 10⁻ (34 h). At this time, the heart tube (HT) is about 500 μm long and has a diameter of about 200 μm (Fig. 1A). It is composed of three layers (Fig. 2A): an outer, two-cell-thick layer of myocardium (MY), a thick middle layer of extracellular matrix forming the cardiac jelly (CJ), and an inner, one-cell-thick endocardium (EN), which surrounds a relatively small lumen. The dorsal side of the HT is open and attached to the body of the embryo through a structure called the dorsal mesocardium (DM), while a membrane called the splanchnopleure (SPL) covers the early HT and contacts its ventral surface (Fig. 2A). In addition, the HT is connected to a pair of omphalomesenteric veins (OVs, sometimes referred to as “veins” in this paper) at its caudal (lower) end and to the conotruncus (outflow tract) at its cranial (upper) end (Fig. 1). The first myofibrillar contractions occur soon after HT formation,8 but effective blood flow does not begin until stage 12 (48 h).20,22

FIGURE 2.

(A) Schematic of cross-section of stage 10⁻ heart tube in the chick embryo (from Manasek14). MY= myocardium; CJ=cardiac jelly; DM=dorsal mesocardium; SPL=splanchnopleure; EN=endocardium. (B) Cross-section of heart tube used in the computational model.

During c-looping, the HT deforms into a c-shaped tube by stage 12. During this process, ventral bending and dextral torsion send the original ventral side of the straight HT to the outer curvature of the looped HT, while the original dorsal side becomes the inner curvature (Fig. 1). After stage 11 (42 h), the dorsal side of the HT gradually closes and the DM ruptures, leaving the tube attached to the embryo only at its ends. In addition, the OVs gradually fuse to add new segments, the primitive left ventricle and atrium, to the caudal end of the HT.

COMPUTATIONAL METHODS

A computational model was constructed for the stage 10⁻ embryonic chick heart. The model captures morphogenesis during the crucial initial phase of c-looping, when left-right directionality is established. During this period (from stage 10⁻ to stage 10⁺, or about 34–38 h of incubation), torsion is not yet pronounced, and the SPL likely does not contribute significantly to the deformation. This membrane, therefore, was ignored in the present model.

The model was created using the commercial program ABAQUS (ABAQUS, INC). A large deformation analysis was used with material properties and internal morphogenetic forces included through a user subroutine (see below).

Model Geometry

The computer-aided design package SolidEdge (EDS software) was used to construct the model geometry, which is based on microscopic images like those shown in Figs. 1 and 5A. The HT and OVs are modeled, respectively, as a straight tube and a pair of curved tubes. The following features of a typical stage 10⁻ heart are included (Fig. 3A,B): (1) The HT has a tapered circular cross section with the smallest diameter located at the conotruncus and largest at a point slightly above the junction of the HT and OVs (Figs. 1A and 5A); (2) the average diameter of the OVs is roughly half that of the HT (Fig. 1A); and (3) the diameter of the left OV is slightly larger than that of the right OV (Fig. 1A,B). The conotruncus is modeled as an extension of the HT.

FIGURE 3.

Computational model for stage 10⁻ chick heart. (A) Front (ventral) view showing regional morphogenetic mechanisms. Also shown are the local circumferential (2) and longitudinal (3) directions for the OVs (dotted lines). HT=heart tube, OV=omphalomesenteric vein, C=conotruncus. (B) Side view showing circumferential (2) and longitudinal (3) directions for the HT (dotted lines). Regions with fixed boundary conditions are indicated in (A) and (B). (C) Circumferential distribution of stretch ratios for actin polymerization in the HT.

FIGURE 4.

Comparison of experimental and model-predicted ventral bending for the heart tube (HT). Side view (ventral side toward the left) of isolated HT and model at stage 10⁻ (A,A′) and stage 10⁺ (B,B′). During the initial stages of looping, the ventral surface (left side of image) bulges outward, while the dorsal surface (right side of image) is relatively straight (B,B′). In (A) and (B), the dark tissue along the dorsal side of the HT is the dorsal mesocardium. Note: Hearts in (A) and (B) are from different embryos.

FIGURE 5.

Comparison of experimental and model-predicted torsion during early looping (ventral view). Chick heart and model are shown at stage 10⁻ (A,A′) and stage 10⁺ (B,B′). During normal c-looping, labels placed near the midline of the stage 10⁻ heart (A) migrate rightward toward the outer curvature at stage 10⁺ (B). Midline nodes in the model (A′, marked by black circles) also migrate rightward (B′). The distances L ₁, L ₂, l ₁, and l ₂ were used to compute the rotation coefficient (see text); dotted lines in (A,B) bound the region where the rotation coefficient was measured. Asterisks in (A,B) denote the position of a stationary label placed near the first pair of somites.

The cross sections of the HT, conotruncus, and OVs are similar. In the model, all three consist of an outer myocardial layer (presumptive myocardium in the OVs) and an inner layer of CJ (Fig. 2B). The relatively thin endocardium is not included. Based on sectional histology data (not shown), the thicknesses of the MY and CJ are taken as 18 and 56 μm, respectively. For boundary conditions, the far ends of the OVs and the dorsal sides of the DM and conotruncus are fixed (Fig. 3A,B).

The created geometry was imported into the meshing program, Patran (MSC Software), as ParaSolid Transmit files and discretized using second order tetrahedral elements. The final mesh contained about 60,000 elements and 198,000 nodes. (The exterior mesh can be seen in Fig. 5A′,B′; the cross section is shown in Fig. 2B.)

Material Properties

Material properties are based on the measurements for the stage 12 chick heart by Zamir and Taber.33 Unfortunately, no data are yet available for earlier stages. These authors determined pseudoelastic strain-energy density functions for the MY and CJ.

The CJ is a passive material, but the MY can undergo active contraction. Because heartbeat is not required for looping to occur,11,17 only cytoskeletal contraction is included in our model. For both MY and CJ, the passive strain-energy density function was taken in the form

$$W_{\rm P} = \frac{A}{B}\left( {{\rm e}^{B(I_1 - 3)} - 1} \right)$$

(1)

where A and B are material coefficients, and $I_1$ is the first strain invariant.

Zamir and Taber33 did not separate passive and active behavior in their analysis, but Zamir et al.32 estimated approximate material parameters for transversely isotropic myocardium. With $W_{\rm P}$ as given above, the active strain-energy density function is given by

$$W_A = C_a (t)\big(\lambda _{\rm f}^2 - 1\big)^2$$

(2)

where $\lambda _{\rm f}$ is the fiber stretch ratio (in the contraction direction) relative to the active zero-stress state (see below). For cytoskeletal contraction, the active material constant was taken as an exponential function of time, $C_a (t) = C_{{\rm max}} (1 - {\rm e}^{ - at} )$, where a is a constant.

As a first approximation, the form of the passive strain-energy function for the OVs was taken to be the same as that of the HT. Also, the properties of the DM and conotruncus were set to be the same as those of the CJ and HT, respectively. The passive parameter values used in the model are the following:

$$A_{{\rm MY}} = 11.1\,{\rm Pa}\quad B_{{\rm MY}} = 0.49\\ A_{{\rm CJ}} = 3.0\,{\rm Pa}\quad B_{{\rm CJ}} = 0.10.$$

(3)

In regions where cytoskeletal contraction was specified, we took

$$C_{{\rm max}} = 50\,{\rm Pa}{\rm .}$$

(4)

In non-contracting regions, we set $C_a (t) = 0$ for all t. We found that the simulation results were relatively insensitive to reasonable changes in material properties. Local material orientations were defined through the user subroutine ORIENT in ABAQUS. Orientation for the HT was defined using a tapered cylindrical coordinate system, while a toroidal coordinate system was used for the OVs. Circumferential and longitudinal directions for the HT and OVs are shown in Fig. 3A,B.

Modeling of Morphogenetic Processes During Looping

The mechanics of many morphogenetic processes can be modeled effectively using tissue growth and contraction algorithms. The theory is well described in the literature,21,26,27 and only a brief summary is presented here. Growth, actin polymerization (cell-shape change), and active contraction are modeled by changing the zero-stress configuration for each material element through the deformation gradient tensors ${\bf F}_g ,{\bf F}_p$, and ${\bf F}_a$, respectively. Cell-shape changes and contraction are assumed to occur with no change in volume; hence det ${\bf F}_p = {\rm det}\,{\bf F}_a = 1$, but ${\rm det}\,{\bf F}_g \ne 1$ in general. The passive and active Cauchy stress tensors are given by

$${\sigma} _{\rm P} = ({\rm det}\,{\bf F}_{\rm P} )^{ - 1} {\bf F}_{\rm P} \cdot \frac{{\partial W_{\rm P} }}{{\partial {\bf E}_{\rm P} }} \cdot {\bf F}_{\rm P}^{\rm T}\\ {\sigma} _{\rm A} = ({\rm det}\,{\bf F}_{\rm A} )^{ - 1} {\bf F}_{\rm A} \cdot \frac{{\partial W_{\rm A} }}{{\partial {\bf E}_{\rm A} }} \cdot {\bf F}_{\rm A}^{\rm T} ,$$

(5)

where T denotes the transpose. With the total deformation gradient tensor relative to the original reference configuration given by F, the tensors ${\bf F}_{\rm P} = {\bf F} \cdot {\bf F}_g^{ - 1} \cdot {\bf F}_p^{ - 1}$ and ${\rm F}_{\rm A} = {\bf F} \cdot {\bf F}_g^{ - 1} \cdot {\bf F}_p^{ - 1} \cdot {\bf F}_a^{ - 1}$ describe the deformation relative to the passive and active zero-stress states, respectively.26

The key to modeling cardiac morphogenesis lies in the tensors ${\bf F}_g ,{\bf F}_p$, and ${\bf F}_a$. Let the unit vectors ${\bf e}_\theta$ and ${\bf e}_z$ be directed along the principal material directions (e.g., circumferential and longitudinal fiber directions), parallel locally to the surface of the heart, with ${\bf e}_r$ being normal to ${\bf e}_\theta$ and ${\bf e}_z$. For a first approximation, we ignore shear components and take

$${\bf F}_g = \lambda _{gr} {\bf e}_r {\bf e}_r + \lambda _{g\theta } {\bf e}_\theta {\bf e}_\theta + \lambda _{gz} {\bf e}_z {\bf e}_z$$

$${\bf F}_p = \lambda _{pr} {\bf e}_r {\bf e}_r + \lambda _{p\theta } {\bf e}_\theta {\bf e}_\theta + \lambda _{pz} {\bf e}_z {\bf e}_z$$

$${\bf F}_a = \lambda _{ar} {\bf e}_r {\bf e}_r + \lambda _{a\theta } {\bf e}_\theta {\bf e}_\theta + \lambda _{az} {\bf e}_z {\bf e}_z$$

(6)

where the $\lambda {\rm s}$ are stretch ratios. Our models simulate morphogenetic mechanisms as follows:

• 1. Growth: Change in the $\lambda _{gi}$ without change in mechanical properties.

• 2. Cytoskeletal contraction: $\lambda _{ai} < 1$ along principal actin fiber directions, with constant element volume $(\lambda _{ar} \lambda _{a\theta } \lambda _{az} = 1)$ and increasing active modulus $C_a$.

• 3. Actin polymerization: $\lambda _{pi} > 1$ along principal actin fiber directions, with constant element volume $(\lambda _{pr} \lambda _{p\theta } \lambda _{pz} = 1)$ and no change in properties.

• 4. Cell migration: The migration of cells from region A to region B is modeled (to a first approximation) as simultaneous negative growth in region A and positive growth in region B.

These processes were included in the model via the ABAQUS subroutine UMAT, which has access to the total deformation gradient tensor F. With ${\bf F}_g ,{\bf F}_p$, and ${\bf F}_a$ specified, the tensors ${\bf F}_{\rm P}$ and ${\bf F}_{\rm A}$ were computed as above and used in the constitutive relations (5). Then, the total stress is ${\sigma} = {\sigma}_{\rm P} + {\sigma}_{\rm A}$.

The free variables in the model are the growth, polymerization, and activation stretch ratios, $\lambda _{gi} ,\lambda _{pi}$, and $\lambda _{ai}$. These parameters were specified regionally according to our hypothesis for c-looping (see below), and their values were adjusted by trial and error until the experimental and model-predicted global shape changes and local morphogenetic stresses and strains matched.

Morphogenetic Force Distributions

In the model, the DM and conotruncus are modeled as passive tissues that neither grow nor contract. In other regions, internal loads were specified according to the experimental observations described below (see Fig. 3A).

Cell Shape Changes Mediated by Actin Polymerization Drive the Bending Component of Looping

Experimental perturbations of actin13 and observed changes in actin orientation and myocardial cell shape during looping10,15,23 suggest that circumferential and circumferential-longitudinal actin polymerization occurs in the dorsal and ventral parts of the myocardium, respectively (see Discussion Section). Consistent with these data, values of the $\lambda _{pi}$ in the HT were specified to vary smoothly from the dorsal side to the ventral side (Fig. 3C). Here, stretch ratios greater than unity denote cell elongation due to actin polymerization, and stretch ratios less than unity denote cell shortening to preserve volume. Note that only circumferential lengthening was specified along the original dorsal midline, and equal circumferential and longitudinal lengthening was specified along the original ventral midline (Fig. 3C). The shape of the radial polymerization curve is a consequence of the constant volume constraint.

The Omphalomesenteric Veins Play a Prominent Role in Initiating Dextral Torsion

In the OVs, progressive longitudinal contraction and growth were stipulated along the caudal and cranial regions, respectively. The contraction is based on data from our laboratory indicating that tension on the caudal side of the veins relaxes when the OVs are exposed to cytochalasin B.30 Growth on the cranial side represents migration of cells through the veins toward the HT.4

Cardiac Jelly Volume Increases During c-Looping

Cardiac jelly was assumed to grow isotropically to simulate experimentally observed accretion and swelling of this extracellular matrix.16,19

The looping simulation begins at stage 10⁻ with $\lambda _{gi} = \lambda _{pi} , = \lambda _{ai} = 1$ everywhere for $i = 1,2,3$. The values determined for these stretch ratios at stage 10⁺ are given in Table 1; the values at stage 10 are midway between 1 and the corresponding value in the table. Components not specified in the table were set to unity. For instance, the HT and the cranial half of the OVs do not contract and, therefore, $\lambda _{ar} = \lambda _{a\theta } = \lambda _{az} = 1.0$ in these regions.

TABLE 1. Growth, actin polymerization, and active contraction stretch ratios used in the model.

The intrinsic loads were applied in discrete steps. Typically, a growth stretch ratio of 1.5 was applied in five to six time steps. All simulations were run on a personal computer with a Pentium IV processor with 1 GB of RAM. Once a solution was obtained, results were analyzed using ABAQUS/Viewer.

EXPERIMENTAL METHODS

Embryo Culture

Fertile white Leghorn chicken eggs were incubated in a forced draft incubator at 38.5°C and constant humidity. We studied c-looping during stages 10⁻ to 11 (approximately 34–42 h of a 21-day incubation period). Embryos were cultured using our recently developed method that eliminates surface tension by submerging the embryos under a thin layer of fluid.31 Briefly, a filter paper ring (Whatman # 2) was used to remove the embryo and vitelline membrane from the egg. The embryo was then sandwiched using a second filter paper ring and placed ventral side up in a Petri dish containing 89% Dulbecco’s modified Eagle’s medium (DME), 10% chicken serum, and 1% antibiotics. To keep the embryo submerged, a stainless steel ring was placed on the assembly. The Petri dish was placed in a plastic bag filled with a mixture of 95% O₂ and 5% CO₂, and the sealed bag was put into a culture incubator, where the embryo develops normally under control conditions.31

Local Longitudinal Strain Measurements

Morphogenetic strains during looping were measured in specified regions of the myocardial layer of the HT and OVs by tracking the motions of fluorescent tissue labels.

Tissue Labeling

Small groups of cells in the myocardium were labeled using a 1:10 dilution of the fluorescent dye Dil (Molecular Probes, D-282) in dimethyl sulfoxide. The dye was injected using a pneumatic pump (World Precision Instruments, PV830) via a glass micropipette controlled by a micromanipulator. Labels were placed about 50 μm apart in linear arrays in different regions (see Fig. 5A,B). Then, the embryo was cultured as described above.

Two types of strains occur in the developing embryonic heart: morphogenetic strains due to developmental shape changes and functional strains due to the heartbeat. Since the heartbeat does not affect looping,11,17 functional strains were eliminated by temporarily stopping the heart by cooling the culture dish and embryo in an ice bath. Studies in our laboratory have shown that the heart stops beating within several minutes of cooling, and development proceeds normally once the embryo is re-warmed.

Intermittently, therefore, embryos were removed from the incubator, staged according to the system of Hamburger and Hamilton,5 and cooled to arrest the heartbeat. This study focused on stages 10⁻, 10, 10⁺, 11⁻ and 11 (approximately 34, 36, 38, 40, and 42 h, respectively). Next, a stack of high magnification sectional images was taken using a fluorescence microscope (Leica DMLB) with a z-spacing of 5 or 10 μm, and the embryo was returned to the incubator. The images were later analyzed for strains (see below), which are reported for seven non-intersecting regions (see Fig. 6).

FIGURE 6.

Experimental and model-predicted longitudinal Lagrangian strains for seven different regions plotted as a function of developmental stage. Cumulative strains referred to stage 10⁻ are shown. In each plot, the solid and dotted lines denote experimental and model strains, respectively. Plots for the positive strain regions (H1, H2, and V2) have the same axis limits, as do the plots for the negative strain regions (H3, V1, V3, and V4). For all cases, model strains are plotted up to stage 10⁺. Error bars for the experimental data denote SEM.

FIGURE 7.

Experimental and model-predicted average circumferential Lagrangian strains in the heart tube and omphalomesenteric veins (OVs), plotted as a function of developmental stage. In each plot, the solid and dotted lines denote experimental and model strains, respectively. Cumulative strains referred to stage 10⁻ are shown. Diameter measurements in each of the three regions were made in the sections bounded by the dotted lines. Error bars for the experimental data denote SEM.

Great care was taken to ensure that the strain data were as accurate as possible. First of all, spacing between the labels was set at an optimal level. If the labels are too far apart, then the effects of curvature become significant. On the other hand, if the labels are too close together, then optical resolution introduces errors. A detailed analysis of the error is given in Taber et al.29

Second, to access the HT, the splanchnopleure (SPL) had to be cut. For regions H1 and V1 (Fig. 6), a small hole was cut on the right side of the membrane to gain access to the right lateral edge of the heart. Even though the hole was relatively small (see arrows in Fig. 9B,B′) and usually healed quickly, a similar sized hole was introduced on the opposite side to maintain symmetry. Similar remarks apply for regions H3 and V2 along the left lateral edge. For region H2, a small hole was cut in the SPL just above that region of the heart. In regions V3 and V4, the myocardium is bound tightly to the SPL and may be labeled directly. Finally, hearts that deviated significantly from the normal phenotype (e.g., hearts that did not rotate) were rejected.

FIGURE 8.

Experimental and model-predicted stress fields at stage 10⁺. (A) Map of experimental stress field. LT=longitudinal tension, LC=longitudinal compression, CT=circumferential tension. Regions (H1–H3, V1–V4) are reproduced for convenience. (B) Model-predicted longitudinal stress field. (B′) Model-predicted circumferential stress field. Also shown are side views of the heart tube showing longitudinal stress distributions with accretion of cardiac jelly included (C) and excluded (C′). Color scheme for B, B′, C, and C′: blue=compression; gray, red, yellow, green=decreasing amounts of tension. R=right side, L=left side, V=ventral side, D=dorsal side.

FIGURE 9.

Effect of dorsal mesocardium (DM) disruption on torsion in chick heart (ventral view). Fluorescent labels were placed along the left and right lateral edges to visualize rotation. (A,A′) Control heart at stage 10⁻ and 4 h later. Labels move away from the edges as heart rotates. Note that the right lateral labels move behind the heart and are seen through the transparent heart tube. (B,B′) Heart with DM dissected at stage 10⁻ and 4 h later. Labels do not migrate from the edges, indicating absence of rotation. Arrows in B,B′ show typical holes cut in the splanchnopleure for access to the heart.

Strain Analysis

A custom MATLAB program provided three-dimensional coordinates for the selected labels. For a user-defined threshold, the program identifies each label as a volume in the image stack and computes the coordinates of the label centroid. The centroid may lie between the image sections.

The distance between the labels was computed at an initial stage (e.g., stage 10⁻) and at a final stage (e.g., stage 10). The stretch ratio, $\lambda _{10/10^-}$, was computed using these distances and the incremental Lagrangian strain was calculated using

$$E_{10/10^ - } = \frac{1}{2}\big(\lambda^2_{10/10^-}- 1\big)$$

(7)

where $\lambda _{10/10^ - }$ indicates the stretch ratio at stage 10 relative to stage 10⁻. The procedure was then repeated using another pair of labels from the array. When all the labels were exhausted in a given embryo, the strains were averaged, giving one data point for a given region. At least six data points were collected for each interval (stage 10⁻ to 10, 10 to 10⁺, 10⁺ to 11⁻, and 11⁻ to 11). In total, approximately 100 embryos were used to obtain the strain data.

Cumulative strains referred to stage 10⁻ were computed using the incremental stretches. For instance, when $\lambda _{10^ + /10}$ and $\lambda _{10/10^ - }$ were known, the cumulative strain at stage 10⁺, relative to stage 10⁻, was found from

$$E_{10^ + /10^ - } = \frac{1}{2}\left( {(\lambda _{10^ + /10} \lambda _{10/10^ - } )^2 - 1} \right).$$

(8)

Average Circumferential Strain Measurements

Average circumferential strains were measured in the HT, right OV, and left OV. For each embryo (n=7), the diameter of the HT and the two OVs were measured at stages 10⁻, 10, 10⁺, and 11⁻. For each data point, five different measurements were taken and averaged. In the HT, measurements were made roughly in the middle section, away from the conotruncus and the OVs (see Fig. 7). In the OVs, measurements were made in regions fairly close to, but still distinct from, the HT (Fig. 7). All measurements were made on high magnification images using ImageJ (NIH software). The incremental stretch ratio λ=D/D ₀ was computed, where D ₀ and D are the initial and final diameters, respectively. Then, cumulative Lagrangian strains were calculated as described above.

Stress Measurements

Regional myocardial wall stresses were characterized via tissue dissection. Small cuts were made in the HT and OVs using fine glass microneedles (n=7 in the HT and n > 50 in the OVs). Visual inspection of the wound indicated the qualitative nature of the stress field in a given region. After cutting, the wound opened or closed, indicating that the tissue was in tension or compression, respectively, in the direction perpendicular to the cut. Tension, compression, or no significant stress was recorded. Because of the manner in which the SPL partially envelopes the OVs,30 longitudinal cuts could not be made in the OVs, and only longitudinal stresses are reported here. Despite the large number of embryos used, the stresses in every region were relatively consistent across all embryos.

Quantification of Dextral Rotation

A rotation coefficient R was computed to quantitatively compare the amount of rotation given by the model with that from experiments. The definition of R is based on the lateral motions of labels along the ventral midline of the HT relative to the left and right edges of the heart (n=7). Let L ₁ and L ₂ represent the distances of the labels from the left and right edge of the stage 10⁻ HT, respectively, with l ₁ being l ₂ the corresponding distances at stage 10⁺ (see Fig. 5A,B). Then, we define

$$R = \frac{{l_2 /l_1 }}{{L_2 /L_1 }} - 1,$$

(9)

which is a positive (negative) number for rightward (leftward) rotation. Rotation at the top and bottom of the HT is limited by attachments to the conotruncus and the OVs, respectively. Hence only those labels in the middle one third of the HT were used to compute the rotation coefficient (Fig. 5). Similar remarks apply to the model, where nodes on the midline of the HT were used to compute R.

Perturbation Studies

In addition to the HT, the other major components of the stage 10⁻ heart include the conotruncus, the DM, and the two OVs. As all of these structures can affect looping morphology, the effects of each were examined by removing them one at a time. Eliminating the contributions of the OVs and outflow tract was relatively straightforward. The tissues in question were severed from the HT using scissors and glass microneedles (n=5 for conotruncus dissection and n=2 for removal of each of the OVs).

To study the effects of the DM, glass microneedles were inserted at the dorsal side of the stage 10⁻ HT through small holes cut in the SPL (n=4). The DM was cut by gently scraping it and simultaneously lifting the heart. This procedure was repeated every 2 h to reduce the effects of DM regeneration. Labels were placed along the left and right sides of the HT to help visualize rotation (see Fig. 9).

Statistics

Strains are reported as mean±standard error of the mean (SEM). Computation of the SEM for cumulative strains is based on the theory of propagation of error, as described by Beers.3

RESULTS

Global Morphogenetic Shape Changes

Ventral Bending

Ventral bending in isolated hearts is compared with that from a model of the HT alone, i.e., without the OVs, DM, and the conotruncus (Fig. 4). Bending in the model is driven by the simulated distributions of actin polymerization (Fig. 3C). The overall shape change predicted by the model at stage 10⁺ closely matches that of the hearts excised from stage 10⁺ embryos (n=18). In both the experimental and the model HT, the dorsal side is relatively straight while the ventral side bulges outward (Fig. 4B,B′).

Dextral Torsion

Dextral torsion can be visualized from the rightward movement of labels placed along the ventral midline of the HT. Comparing the motions of experimental labels (Fig. 5A,B) and corresponding nodes in the model (Fig. 5A′,B′) shows that rotation at stage 10⁺ is qualitatively similar between the model and the experiment. In addition, the amount of rotation produced by the model was quantitatively compared with that of the experiment. The experimental rotation coefficient R was 0.499±0.182 (mean±SD, n=7), while the value given by the model was 0.594. The model also predicts that very little rotation is possible when the OVs are passive (R=0.103).

During looping, fusion of the OVs adds new segments to the caudal end of the HT.30 This process, possibly coupled with longitudinal cytoskeletal contraction, results in a perceived downward movement of the caudal end of the heart as looping progresses. This can be seen by noting the position of the bottom of the OVs relative to a stationary label (denoted by the asterisk in Fig. 5A,B) placed near the first pair of somites. Although OV fusion is not included in the present model, a similar downward movement is observed in the model (Fig. 5A′,B′) due to contraction in the caudal OVs. Note that this downward shift is not observed when OV contraction is eliminated from the model (compare Figs. 10C and 10C′).

Local Morphogenetic Shape Changes

Local Longitudinal Strains

Longitudinal strains in the myocardium were measured by tracking labels from stage 10⁻ to stage 11 in seven regions of the heart (Fig. 6). Cumulative strains relative to stage 10⁻ were computed from the incremental strain at each stage.

Experimental (solid line) and model-predicted (dotted line) longitudinal strains are compared in Fig. 6. The experimental strain field was highly nonhomogenous, with strain magnitudes (and sometimes signs) differing considerably between regions. In the OVs for instance, the strains were negative in regions V1, V3, and V4 while region V2 exhibited large positive strains, even reaching a magnitude near 1.0 at stage 11.

In the HT, the experimental strains were positive in regions H1 and H2. Although the mean strain at stage 11 for these two regions was roughly the same (0.465 for region H1, and 0.471 for region H2), the trends were somewhat different. In region H2, the strains increased from stage 10⁻ to stage 10⁺ and then remained relatively constant, while a steady increase was observed in region H1. In region H3, the strains initially increased, followed by a steady decline.

It is also worth noting the wide range of strain magnitudes observed at each stage. Isolating just the negative strains at stage 11⁻, magnitudes ranged from −0.048 (region H3) to −0.358 (region V3). The larger negative strains caused the labels to move so close together that, by stage 11, they could not be distinguished by the image processing software. Hence data are not reported in regions V1, V3, and V4 at stage 11. For positive strains at stage 11, the range was from 0.471 (region H2) to 0.900 (region V2).

Most of the model-predicted strains closely follow the experimental strains (Fig. 6), with one exception being the strain at stage 10⁺ in region H3. Here, the computed strains monotonically increase in magnitude, while the experimental strains reversed direction at stage 10. Note also that strain magnitudes in this region were much smaller compared to those in other regions (Fig. 6).

Average Circumferential Strains

Gross diameter changes in the HT and OVs were measured. Experimental and model-predicted average circumferential strains are compared in three regions of the embryonic heart (Fig. 7).

In the HT, the average strain increased steadily from stage 10⁻ to stage 10⁺ and then changed little at stage 11⁻. In the left OV, the value steadily increased from stage 10⁻ to 11⁻. The strains in the right OV were smaller when compared to those of the other two regions, and did not increase much after stage 10. The diameter could not be measured for stage 11⁻ as the rotated HT obstructs much of the right OV at this stage. The model closely follows the experimental trends in all three regions.

Radial Strains

Manasek et al.15 reported that the myocardium of the HT thickens during looping at the inner curvature and thins at the outer curvature. Relatively crude measurements obtained from sectional images in Manasek et al.15 yielded average radial cumulative Lagrangian strains (referred to stage 10⁻) of 0.212 and −0.314 at the inner (dorsal) and outer (ventral) curvature, respectively. These values are for stage 11⁺ (44 h). The respective model-predicted values at stage 10⁺ (38 h) are 0.110 and −0.284.

Strains from Left-Looping Hearts

During normal morphogenesis, leftward torsion is extremely rare.30 During data collection for longitudinal strains, two hearts with labels in the cranial region of the OVs looped to the left. In these hearts, region V1, which normally had negative strains, had positive strains (E=0.156 at stage 10⁺ and E=0.496 at stage 11⁻, both values relative to stage 10). Similarly, region V2, which normally had positive strains, had negative strains (E=−0.047, −0.078, and −0.395 at stages 10⁺, 11⁻, and 11 respectively, all values relative to stage 10).

Morphogenetic Stresses

Experimental stress fields in the myocardium showed a nonhomogenous pattern at stage 10⁺ (Fig. 8A). In the longitudinal direction, tensile stresses were observed throughout the HT and most of the caudal regions of the OVs, while compressive stresses were found in the cranial parts of the OVs. The longitudinal stresses given by the model matched the experimental pattern of tension and compression (Fig. 8B). In the circumferential direction, the myocardium of the HT was in tension (Fig. 8A,B′).

FIGURE 10.

(A,A′,B,B′): Effects of mechanical perturbations on looping morphology (ventral view). (A,A′) Effect of removing right omphalomesenteric vein (OV) at stage 10⁻. (A) Experimental heart at stage 12. (A′) Model-predicted shape at stage 10⁺. Note the increased curvature of the left OV (arrows). (B,B′) Immediate effect of severing the conotruncus at stage 10⁺. Both the experiment (A) and the model (B) show that the heart tube tilts rightward. Note also that the OV shapes are similar between model and experiment. (C–C″′): Effect of morphogenetic force perturbations on looping morphology in the model. In all cases, actin polymerization in the heart tube is turned off, and only contraction in the caudal OVs is perturbed. (C) Baseline model with contraction specified in both veins. (C′) Contraction turned off in both veins; normal downward motion of the OVs is not observed. (C″) Contraction turned off in left OV only; some downward movement is observed (C″′) Contraction turned off in right OV only; the heart tube moves to the left.

According to the model, the entire myocardium of the HT is in a state of longitudinal tension that increases from the ventral to the dorsal side (Fig. 8C). When CJ inflation is not included, however, the ventral side is in compression, while the dorsal side is in tension (Fig. 8C′). Note also that the normal tension in the myocardium causes compression of the cardiac jelly (results not shown).

Perturbation Studies

Mechanical Perturbations

Figure 9 shows the effect of disrupting the dorsal mesocardium (DM) on dextral rotation. For technical convenience, edge labels, rather than midline labels, were used in this part of the study. In control embryos, dextral rotation caused labels from the left and right sides of the HT to appear to move closer together (Fig. 9A,A′). In contrast, when the DM was dissected, the labels remained fixed on their respective edges, indicating a lack of rotation (Fig. 9B,B′). Four embryos were used; two did not rotate and two rotated normally when the DM was cut.

In the computational model, removing the DM leaves the HT completely free along its dorsal extent, as well as at its top face (see Fig. 3B). These are not realistic boundary conditions, as the embryonic heart in the experiments remained fixed to the foregut near the cranial and caudal ends. Therefore, the conotruncus in the model was extended and fixed at the top to simulate the attachment to the developing aortic arches. In the model with DM removed, rotation was severely diminished (results not shown). The value of the rotation coefficient R for the perturbed model is 0.181, which is a 70% decrease from the value obtained with the model under control conditions (R=0.594).

When the right OV was removed, the longitudinal curvature of the left OV was greater than normal in both the experimental and the model hearts (compare Figs. 1C and 10A and Figs. 5B′ and 10A′). Although the experimental hearts were at a later stage of development than the model, the general shapes of the model HT and OV are similar to those from experiment (Fig. 10A,A′). Removal of the left OV resulted in left looping in both model and experiment (results not shown).

Freeing the top end of the HT and the DM resulted in a rightward tilt of the HT as observed in both experiments and the model (Fig. 10B,B′). Note also that the model-predicted shape of the OVs after conotruncus dissection is similar to that of experiment (Fig. 10B,B′).

Morphogenetic Force Perturbations

To study the effects of OV forces on dextral torsion, simulations were run with morphogenetic forces in the OVs turned off one region at a time. To avoid complications due to the coupling between bending and rotation, ventral bending was eliminated by turning off actin polymerization in the heart tube for all simulations in this section. This led to a rotation coefficient R=0.30 for the baseline case with the only loads being cytoskeletal contraction and cell migration specified in both OVs (Fig. 10C).

When contraction is turned off in both OVs, there is a sharp decrease in rotation (R=0.098 Fig. 10C′). When contraction is turned off in only the left OV, the unbalanced contraction in the right OV causes a sharp increase in rotation (R=0.605, Fig. 10C″). When contraction is turned off in only the right OV, the unbalanced force in the left OV causes slight leftward torsion (R=−0.1, Fig. 10C″′).

It is interesting to note the shape of the OVs predicted by the model under various perturbations. When contraction is eliminated from both OVs, they do not move downward (Fig. 10C′) as in the baseline case (Fig. 10C). With contraction eliminated from the left OV, some downward movement occurs (Fig. 10C″). With contraction eliminated from the right OV, the right OV elongates and the HT shifts leftward (Fig. 10C″′).

DISCUSSION

Over the years, a number of ideas have been put forward to explain the mechanisms that cause bending and torsion during cardiac c-looping.28 Physical and mathematical models have been used to demonstrate the physical plausibility of some of the proposed mechanisms,9,16,19,20 but most of these hypotheses have since been contradicted by new experiments. Recently, we have proposed a new hypothesis for the initial stages of c-looping.13,28,30 According to this theory, changes in cell shape driven by actin polymerization in the HT are responsible for the bending component of c-looping, while unbalanced forces in the OVs, due to a combination of cell migration and cytoskeletal contraction, initiate dextral rotation (torsion). Forces exerted by the splanchnopleure (SPL) cause further torsion. In other words, ventral bending is intrinsic to the HT, while torsion is caused by external forces.

In prior work, this hypothesis has been tested experimentally using chemical and mechanical perturbations. We have found that disrupting actin polymerization using the drugs latrunculin A and cytochalasin D inhibits the bending component of c-looping.13 Moreover, the HT untwists when the OVs are severed and the SPL removed.30 These results support our hypothesis.

In addition, we have previously examined the feasibility of our hypothesis using finite element models.13,30,32 These models feature relatively simple geometries and loads, and most comparisons to experimental results were qualitative in nature. The current study tests the hypothesis quantitatively using a more realistic model and new experimental data on morphogenetic stress and strain distributions. New features of the present model include more realistic geometry, anisotropic material properties, separate passive and active myocardial stresses, and CJ growth.

With realistic morphogenetic force distributions specified according to available experimental data, model-predicted global shape changes (Figs. 5 and 4), local morphogenetic strains (Figs. 6 and 7), and local morphogenetic stresses (Fig. 8) agree relatively well with experimental data. When tested under various perturbed conditions, the model also yields results that agree with experimental results (Figs. 9 and 10). Because of the agreement with a fairly wide range of data, we conclude that the proposed model supports our hypothesis for early c-looping.

Remarks on Ventral Bending

The morphogenetic forces simulated in the HT are based on the work of Latacha et al.13 who showed that actin polymerization likely provides the forces that drive cardiac bending. For the assumed distributions of polymerization shown in Fig. 3C, our model for the HT exhibits global bending geometry that is remarkably similar to that of the experimental HT (Fig. 4). The relatively straight dorsal side at stage 10⁺ is a notable feature of both. Moreover, computed strain and stress fields in the HT agree with our data in most regions (Figs. 6, 7, and 8). Taken together with results from previous studies,28 the present results lend further support to the actin polymerization hypothesis.

The distributions of actin polymerization specified in the HT are based on the actin microstructure observed in the myocardium of the looping chick heart.10,23 During c-looping, actin has a predominant circumferential orientation near the inner curvature (IC) of the HT,10,23 while there is no preferred orientation at the outer curvature (OC). According to our model, actin polymerization causes cells near the IC to elongate circumferentially and shorten in the longitudinal direction (to preserve volume), while cells at the OC elongate in both the circumferential and longitudinal directions. Hence, cells near the OC “push” longitudinally, while those near the IC “pull”. In this way, cells throughout the HT contribute to the bending. Latacha et al. 13 showed that actin polymerization forces on both the dorsal and ventral sides of the HT contribute roughly equally to the bending.

Myocardial strain and stress distributions are caused by the combined effects of bending, torsion, and CJ pressure. Ventral bending alone via the proposed mechanism would produce positive longitudinal strain at the OC (region H2) and relatively small strain at the sides (regions H1 and H3) (Fig. 6). Accretion and swelling of the CJ alone would induce positive strains throughout the HT. The combination of bending and CJ growth yields large positive strains in region H2 and relatively small positive strains in regions H1 and H3 (Fig. 6). Interestingly, the experiments gave moderate positive longitudinal strains in region H1 and strains that become negative in region H3. These strains may be caused by the relatively small amount of rightward bending after stage 10,30 which was not included in the model.

Cardiac jelly pressure has previously been suggested as the driving force for c-looping.16,19 Although this idea has been contradicted by experiments with hyaluronidase,2 the CJ does increase in volume during the early stages of looping.16,19 This growth alters myocardial stresses. Without CJ pressure, active cell-shape changes would produce longitudinal compression at the OC and tension at the IC (Fig. 8C′). Growth of the CJ causes a net added tensile stress in both the longitudinal and circumferential directions throughout the HT (Fig. 8C) that likely enhances the structural integrity of the myocardium.16,19

Finally, we note that the model-predicted stress in the HT at stage 10⁺ is about 10 Pa in both the longitudinal and circumferential directions. This value is consistent with the result of Zamir and Taber,33 who estimated an isotropic stress of 48 Pa at stage 12. Continued growth of the jelly would steadily increase the stress from stage 10⁺ to stage 12.

Remarks on Dextral Rotation

According to Voronov et al. 30 the OVs play a central role in initiating rotation and in determining looping directionality. In the model, the combined actions of cytoskeletal contraction in the caudal regions of the OVs and cell migration into the cranial regions produce global shape changes and rotation, as well as OV stress and strain distributions, that are in good agreement with our measurements (Figs. 5, 6, 7, and 8). Also, the model-predicted torsion when the OVs are passive is significantly smaller compared to that of the unperturbed model (Section “Dextral Torsion”). This is consistent with the results of Voronov et al.,30 who show that the HT has no intrinsic ability to rotate during early c-looping.

According to our hypothesis,28 cell migration in each of the OVs pushes the heart tube toward the opposite side while contraction pulls the HT toward the same side. Accordingly, turning off contraction in just the right OV induces leftward torsion while eliminating contraction in the left OV increases the amount of rightward torsion (Fig. 10C″). The model also predicts that the magnitude of the HT rotation depends on the coupled effects of both cell contraction and cell migration in the OVs (see Fig. 10C–C′′′ and Section “Morphogenetic Force Perturbations”). Turning off either one results in a sharp drop in the rotation coefficient.

Moreover, the direction of rotation depends on the asymmetry in vein geometry. During normal rightward looping, the left OV, which is slightly larger (Fig. 1B,C), dominates and pushes the HT to the right.30 The initial OV asymmetry may be genetic in origin.

The strain distributions in the OVs reflect this complex behavior. Since, the ends of the OVs are fixed to the body of the embryo, the contraction produces negative longitudinal strains (Fig. 6) and tensile longitudinal stresses (Fig. 8) on the caudal side of the OVs (regions V3 and V4). The situation on the cranial side of the OVs is slightly more complicated. Migration of precardiac cells into the HT through the OVs normally would result in positive longitudinal strains and compressive stresses in this region. This pattern is indeed true for the left OV (Figs. 6 and 8; region V2). The right OV, however, is pushed rightward by the larger left vein, causing negative strains and compressive stresses (Figs. 6 and 8; region V1). The longitudinal stress and strain patterns in the OVs at stage 10⁺ are summarized in Fig. 11. Note that the model is able to correctly predict these complex distributions (Figs. 6 and 8). In addition the model correctly predicts that, although both the OVs increase in diameter from stage 10⁻ to 10⁺, the left OV grows at a much faster rate (Fig. 7).

Interestingly, we found that the longitudinal strain pattern during left looping, when the right OV dominates, was opposite of that for normal rightward looping, i.e., positive strains were observed in the right cranial OV (region V1), while negative strains were observed in the left cranial OV (region V2) (see Section “Strains from Left-Looping Hearts“). We also observed a similar trend reversal in the circumferential strains (data not shown). These findings further support the idea that the dominant OV determines looping directionality.

FIGURE 11.

Summary of measured longitudinal stress and strain distributions in the omphalomesenteric veins (OVs) during normal c-looping. Items breaking left-right symmetry are shown enclosed in ovals and interchange positions during abnormal left looping.

Remarks on Strain Measurements

To our knowledge, only two previous studies have reported morphogenetic wall strains during looping. Using data originally presented by Manasek et al.,15 Lacktis and Manasek12 analyzed the movement of iron oxide particles from stage 10 to stage 11⁻ in the embryonic chick heart. These particles were placed roughly in the middle of the HT. However, only one embryo was used for the analysis, and the analysis ignored distances in the z-direction (depth).

Stalsberg and DeHaan25 also used iron oxide particles and measured gross changes in length along the right and left lateral edges of the HT from stage 10 to stage 14. With depth information again ignored, longitudinal stretch ratios along the right and left lateral edges of the HT are reported.

In this study, we measured longitudinal and circumferential strains from stage 10⁻ to 11. Since bending and torsion occur simultaneously in the looping heart, the strain field is highly inhomogeneous. Hence, we report strains for seven non-intersecting regions (Fig. 6). Depth information was taken into account by taking a stack of high magnification sectional images. Also, our use of fluorescent labeling overcomes a limitation in using iron oxide particles, as they could be injected without removing the SPL, which influences looping.30,31

The magnitudes of HT strains reported in the present study are in general agreement with those measured by Lacktis and Manasek.12 The overall changes in length are consistent with the data in Stalsberg and DeHaan.25 Also, sectional images in Manasek et al. 15 indicate that the outer curvature of the HT shows a prominent reduction in thickness from stage 10⁻ to stage 11⁺, while the inner curvature shows a modest increase in thickness. The model correctly predicts this trend (see Section “Radial Strains”).

Limitations

In the present model, the algorithms for modeling actin polymerization and cell migration are relatively simplistic, and viscoelastic effects are ignored. Moreover, the current model does not include a mechanism for dextral bending. Future studies should address these issues.

The results from our experiments involving DM dissection also warrant further study. Those experiments yielded an all-or-nothing response with respect to dextral torsion; i.e., either the heart did not rotate or rotation appeared to be normal. In contrast, the model predicts a sharp drop in rotation when the DM is removed. Dextral torsion is an extremely conserved process during development, and other processes may compensate when mechanical boundary conditions are perturbed. This could be the reason why 50% of our dissected hearts looped normally. Another reason may be incomplete dissection, as the limited visual accessibility of the stage 10⁻ DM made it difficult to confirm the extent of the cut or whether the cut healed during incubation.

Our model does not include the SPL, which supplies much of the force that drives cardiac torsion during late c-looping.30 Between stages 10⁻ and 10⁺, however, the amount of rotation is relatively small, and the SPL likely has little effect. In fact, if the membrane is removed from a stage 10⁻ heart, development proceeds normally until stage 10⁺ (results not shown). The current model also does not include OV fusion, which lengthens the caudal end of the heart tube. To model looping beyond stage 10⁺, both the SPL and OV fusion will need to be added.

Ideally, non-collinear arrays of tissue labels should be used to obtain three components of surface strain.1,7 Unfortunately, technical problems prevented us from making these measurements at this time. In particular, we found it extremely difficult to distribute labels properly without excessively damaging the SPL, which envelopes the heart. Hence, local circumferential and shear strains were not measured. However, average circumferential strains were measured by considering gross diameter changes.

In conclusion, we have used a computational model to show that active cell-shape changes can cause ventral bending during c-looping in the embryonic heart. We have similarly shown that unbalanced forces in the omphalomesenteric veins can initiate dextral rotation. Even though it is possible that different sets of morphogenetic forces could produce similar global shape changes, it is highly unlikely that they would also produce the same local stress and strain distributions and perturbation outcomes. Because agreement between experiment and theory is reasonable on all of these fronts, the results from our model support our recent hypothesis for c-looping.28 The model also provides fundamental new insights into this complex problem.