Tunicates, small invertebrates within the phylum Chordata, possess a robust tubular heart which pumps blood through their open circulatory systems without the use of valves. This heart consists of two major components: the tubular myocardium, a flexible layer of myocardial cells that actively contracts to drive fluid down the length of the tube; and the pericardium, a stiff, outer layer of cells that surrounds the myocardium and creates a fluid-filled space between the myocardium and the pericardium. We investigated the role of the pericardium through in vivo manipulations on tunicate hearts and computational simulations of the myocardium and pericardium using the immersed boundary method. Experimental manipulations reveal that damage to the pericardium results in aneurysm-like bulging of the myocardium and major reductions in the net blood flow and percentage closure of the heart's lumen during contraction. In addition, varying the pericardium-to-myocardium (PM) diameter ratio by increasing damage severity was positively correlated with peak dye flow in the heart. Computational simulations mirror the results of varying the PM ratio experimentally. Reducing the stiffness of the myocardium in the simulations reduced mean blood flow only for simulations without a pericardium. These results indicate that the pericardium has the ability to functionally increase the stiffness of the myocardium and limit myocardial aneurysms. The pericardium's function is likely to enhance flow through the highly resistive circulatory system by acting as a support structure in the absence of connective tissue within the myocardium.
Many animals contract flexible, muscular tubes to drive viscous fluid through smaller, resistive tubes, including hearts in several phyla (Xavier-Neto et al., 2007, 2010), drinking in butterflies and moths (Krenn, 2010), lymphatic pumps (Gashev, 2002), and tracheal systems in terrestrial insects (Greenlee et al., 2013; Harrison et al., 2013). These pumps are versatile and widespread, operating over several orders of magnitude of sizes and flow velocities and occurring commonly and independently throughout bilaterians (Xavier-Neto et al., 2007).
Despite being both common and widespread, the dynamics of driving fluid with a flexible boundary that deforms as a result of viscous resistance are not well understood. In typical mathematical models of pumping, the motion of the pump is prescribed and the resulting fluid flow is determined analytically or numerically. Furthermore, accessory structures to these pumps, such as trabeculae, blood cells and external tissues, are often ignored. Here, we focused on the interaction between the flexible, contractile boundary responsible for driving fluid flow and an accessory structure, the pericardium, through the use of animal experiments and a computational model.
Tunicates (urochordates) are small marine invertebrates belonging to the phylum Chordata (Fig. 1A). Recently, tunicates have been investigated as a model system for vertebrate embryonic heart development because: they represent vertebrates' closest invertebrate relatives and are a possible sister taxon of Vertebrata (Delsuc et al., 2006; Davidson, 2007; Lemaire, 2011); they demonstrate a deep phylogeny in the genetic mechanisms of heart development (Nunzi et al., 1979; Burighel et al., 2001; Xavier-Neto et al., 2010); and fluid flow in their circulatory system is driven by a valveless, tubular heart that resembles the size, morphology and kinematics of the hearts of developing vertebrate embryos (Ruppert and Carle, 1983; Dooley and Zon, 2000; Manner et al., 2010; Lemaire, 2011). Thus, understanding the flow produced by the valveless, tubular hearts of tunicates is key to understanding both the ontogenetic and evolutionary development of closed, high-pressure circulatory systems of vertebrates (Santhanakrishnan and Miller, 2011; McMahon, 2012; McMahon et al., 2011; Dooley and Zon, 2000; Moorman and Christoffels, 2003; Xavier-Neto et al., 2007, 2010; Simoes-Costa et al., 2005).
The circulatory system of most tunicates is characterized as open, with major vessels that lead away from the heart and open into sinuses to deliver blood to body tissues. Unusually robust for a sessile, suspension-feeding animal (Xavier-Neto et al., 2010), the tubular heart drives fluid flow through these resistive tissues of the body and reverses flow periodically to even the transport of oxygen, nutrients and waste throughout its circulatory system (Heron, 1975).
The tunicate heart (Fig. 1B) consists of a layer of myocardial tissue formed by an invagination of the outer, stiff pericardium, each one cell-layer thick (Kalk, 1970). The myocardium and pericardium are connected at the raphe on the basement membrane of the pericardium where the invaginated tube has fused together and runs along the ventral surface of the pericardium (Xavier-Neto et al., 2010). The myocardium consists of a field of myoepithelial contractile cells, which are much longer than they are wide, and a 1–2 μm thick basal membrane that lines the lumen (Kalk, 1970). The fusiform cells do not form fibers, but the myofibrils run along the cell length and are at a diagonal in relation to the basement membrane of the pericardium (Kalk, 1970). The pericardium is innervated early in development from the right subendostylar nerves and latero-posterior mantle nerves through the raphe, which possibly have a role in sensing pressure within the cardiac lumen (Bone and Whitear, 1958; Nunzi et al., 1979; Burighel et al., 2001).
List of symbols and abbreviations
- threshold potential
- diffusion coefficient of electropotential
- duration of applied current
- Young's modulus
- flexural stiffness of myocardium
- force per unit area
- force per unit length
- heart rate
- fstr(r,t) and fbend(r,t)
- forces due to resistance to stretching and bending
- forces due to holding target points in place
- second moment of area
- amplitude of applied current
- immersed boundary method
- applied current representing the action of the pacemaker
- resting length of springs
- maximum effective stiffness of springs
- spring stiffness
- stiffness coefficient proportional to resistance to stretching or compression
- stiffness coefficient
- length of heart tube
- diameter of contracted myocardium
- diameter of relaxed myocardium
- pulse number
- fluid pressure
- pericardium diameter
- pericardium-to-myocardium diameter ratio
- Lagrangian position parameter
- correlation coefficient of linear regression
- period of pulse
- fluid velocity
- dimensionless electropotential
- dimensionless blocking potential
- Womersley number
- position along heart tube
- position in fluid grid
- Cartesian coordinates at time t of the material point labeled by r
- prescribed position of target boundary
- y-intercept of linear regression
- resetting rate
- two-dimensional Dirac delta function
- strength of blocking
- fluid dynamic viscosity
- fluid density
- heart beat frequency
The myocardium contracts to reduce the radius of the heart's lumen. This contraction propagates down the length of the heart tube to drive blood flow through the circulatory system. These contraction waves are produced by nodes on either end of the contractile region of the heart tube (Anderson, 1968; Heron, 1975). The pericardium consists of an outer layer of cells, which are stiff and resistant to manipulation. Cells of the pericardial tissue are square and junctions between cells provide mechanical support, although no additional specialized connective tissue exists in the myocardium or pericardium (Kalk, 1970).
For many vertebrates, the pericardium is an important structure. It serves to support the heart in the thoracic cavities of mammals, to limit ventricular fill, and to prevent rotation during ventricular contraction that is associated with higher mortality (Fritz et al., 2012). In elasmobranch and teleost fish, the pericardium supports sub-ambient pressures in the pericardial cavity that surrounds the atrium and ventricle, which aids in filling the ventricles vis a fronte (Sudak, 1965a,b; Shabetai et al., 1979; Farrell et al., 1988). Recently, fluid motion generated by heart beats between the pericardium and epicardium during vertebrate embryonic development was investigated and found to move important precursor cells that contribute to epicardium morphogenesis in zebrafish (Peralta et al., 2013).
Although an important structure for heart function in vertebrates, little is known about the function of the pericardium in tunicates. Because flows in the tunicate heart are viscous dominated, the fluid dynamics of transport are very different from those in the large hearts of most adult vertebrates but are closer to those of embryonic vertebrate hearts. Anderson (1968) and Kalk (1970) posited that the pericardium likely plays an important mechanical role in blood flow after observing that a damaged pericardium resulted in distention of the myocardium and an interruption in contractions. Xavier-Neto et al. (2010) have suggested that distention of the contractile layer is a distinct problem for valveless tubular hearts in general. The pericardium could improve some of the mechanical issues surrounding valveless, tubular hearts by physically restricting distention of the tube.
To investigate the function of the pericardium in tunicates, computational modeling can be a valuable complement to experimental manipulation. Models can be used to eliminate confounding biological features of flow, e.g. chemical signaling, as well as decouple physical features of flow that may contribute differently to biological function. For example, Taber et al. (2007) used a computational fluid dynamics model to show that cardiac cushions helped to initiate the transition between peristaltic and pulsatile flow in tubular hearts. Miller (2011) and Santhanakrishnan et al. (2009) used numerical simulations of fluid–structure interactions to quantify the flow and vortex formation in the early embryonic ventricle. Computational models are also being used to investigate the pumping mechanism by which valveless tubular hearts drive flow (Baird et al., 2014; Loumes et al., 2008; Avrahami and Gharib, 2008).
Mathematical models of pumping can be divided into those in which (1) the motion of the boundary is prescribed and the fluid is driven by that motion, and (2) the fully coupled fluid–structure interaction problem of a pump that drives the flow and is also deformed by the flow is solved (Taber et al., 2007; Jung and Peskin, 2001; Lee et al., 2009, 2012; Loumes et al., 2008; Fritz et al., 2012; Avrahami and Gharib, 2008). The immersed boundary method (IBM) provides a relatively easy way to solve the fully coupled fluid–structure interaction problem (Peskin, 2002) and has been used by several groups to simulate flow in tubular hearts (Baird et al., 2014; Jung and Peskin, 2001; Lee et al., 2009; Jung et al., 2008; Lim and Jung, 2010). Despite its prominence as an anatomical structure, none of these numerical studies have included a rigid structure surrounding the contracting region of the tubular heart.
In this study, we investigated the role of the pericardium in producing fluid flow in the valveless, tubular heart of the solitary tunicate Ciona savignyi Herdman. We used two approaches: (1) experimental manipulation: measurements of blood flow in the lumen and contraction kinematics were taken on animals with in vivo manipulations of the pericardium and otherwise intact circulatory systems, and (2) computational modeling: IBM models were used to simulate flow through a racetrack circulatory system driven by a valveless, tubular contracting region with and without a rigid pericardium.
Experimental manipulation of animals
Relaxed body lengths of tunicates (N=16 animals) varied from 25 to 55 mm. Neither myocardium diameter (F1,15=2.62, corrected P=1) nor pericardium diameter (F1,15=4.75, corrected P=0.50), a measurement that should not change with damage to the pericardium, had a significant relationship with body length. Pericardium-to-myocardium diameter (PM) ratio (see Materials and methods for full definition) also did not have a significant relationship with body length (F1,15=1.47, corrected P=1).
Experimental observations of myocardium kinematics showed that the myocardium becomes inflated within the pericardial space after the pericardium is damaged and pericardial fluid is allowed to leak out of the pericardial space. The average PM ratio dropped from 2.4±0.09 for the undamaged condition to 1.8±0.2 for the damaged condition, and this decrease in PM ratio with increasing damage level was significant (F1,15=80.1, corrected P=2×10−6). We observed that severe damage to the pericardium often resulted in kinks to the myocardium, occluding flow in the lumen during pumping.
After pericardium damage, it was more difficult for the myocardium to contract, but the rate at which the myocardium contracted did not change. The contraction ratio (diameter of the myocardium during contraction to diameter of the relaxed myocardium) ranged from 0.135 (86% closure) to only 0.566 (43% closure) (Fig. 2C). This decrease was not significantly correlated with pericardium damage (F1,14=9.58, corrected P=0.10), but it was significantly correlated with PM ratio (F1,14=12.7, corrected P=0.034). Heart beat frequency (fH) did not change significantly with damage (Fig. 2A; F1,14=0.112, corrected P=1) or PM ratio (F1,14=0.2922, corrected P=1). Actuation of heart beats was not affected by changes in myocardium diameter, but the ability of the myocardium to contract fully was affected.
Peak dye speed in the lumen was severely interrupted by damage to the pericardium (Fig. 3, Fig. 2B,D). There was a proportional response to damage; a significant linear relationship between PM ratio and peak dye speed is shown in Fig. 2D (F1,13=78.3, P=8×10−6). Contraction ratio also negatively correlated with peak dye speed (F1,9=5.27, P=0.002).
Two outliers were identified and omitted. The first was in the peak dye speed measurement (damaged group, PM ratio=1.63, peak dye speed=7.13 mm s−1; Bonferroni P=0.012). We observed that during the dissection, the myocardium near the terminus of the heart (at the point of injection) was inadvertently damaged, which allowed a large amount of blood to leak out. This apparently resulted in an increase in peak dye speed likely due to relief of pressure in the lumen rather than being a result of active pumping. The second outlier omitted was in the contraction ratio measurement (damaged group, PM ratio=1.96, contraction ratio=0.31; Bonferroni P=0.004).
A numerical model of an elastic heart tube connected to a nearly rigid racetrack was constructed to study flow through the Ciona heart. The racetrack design was used for easy comparison to previous models of tubular heart pumping (Jung and Peskin, 2001; Hickerson et al., 2005; Baird et al., 2014; Avrahami and Gharib, 2008; Lee et al., 2012). The contracting section of the tube resists bending and stretching. The racetrack section was constructed by connecting two sections of straight tube (one of which represents the elastic heart) to curved sections. The resting diameter of the racetrack was constant throughout its length. Dimensions and elastic properties of the racetrack are given in Table 1. A diagram of the model set up is shown in Fig. 4.
The contraction waves were driven by a simplified model of the heart electromechanics. The propagation of the action potentials along the heart tube were modeled using the FitzHugh–Nagumo equations (FitzHugh, 1955; Nagumo et al., 1962). The FitzHugh–Nagumo equations are commonly used to model problems in electrophysiology, and they were derived by reducing the Hodgkin–Huxley equations to two variables. Assume that v(x,t) is the dimensionless electropotential that is a function of position along the tube, x, and time, t, such that v=0 is the resting potential and v=1 is the maximum electropotential. Also assume that w represents a dimensionless ‘blocking mechanism’ that represents the action of the potassium channels and the blocking of the sodium ion channels. w=0 is the resting strength, and w=1 is the maximum strength. The 1D FitzHugh–Nagumo equations may be written as: (1) (2)
where D is diffusion coefficient of the electropotential, a is the threshold potential, γ is a resetting rate, ε is the strength of blocking and I(x,t) is an applied current that represents the action of the pacemaker. The equation for the applied current is given as: (3)
where Ltube is the length of the heart tube, Imax is the amplitude of the applied current, n is an integer representing the pulse number, T is the period of the pulse and dp is the duration of the applied current.
The next step is to connect the electrophysiology model with the generation of contraction force. We used a simple model to drive the motion of the 2D heart with 1D muscles. Linear springs with variable spring stiffness and zero resting lengths were used to connect the top and bottom of the tube. The resting lengths of these springs were given as functions of the local electropotential according to the following rule: (4)
where kmax is the maximum effective stiffness of the springs as 0<v<1, and no force is applied by the muscle springs when v<0. kmax was set to ks/100, where ks is spring stiffness. This is only a toy model of muscle activation and force generation, but it has the necessary effect of creating a contraction wave that follows the depolarization along the heart tube.
For each simulation time step, the fluid velocity and pressure profiles were recorded at the outflow end of the myocardium at points across the diameter of the rigid circulatory tube. Two flow speeds were then calculated from the fluid velocity profiles: peak flow speed and mean flow speed. For each time step, the maximum speed obtained in each fluid velocity profile was recorded. Of these, the local maximums for each heart beat were chosen and then time averaged to find the peak flow speed. Mean flow speed was found by spatially averaging the speed across the tube's diameter and then taking the time-averaged mean of those speeds.
From the pressure profiles, four values of pressure were calculated. Pressure across the diameter of the racetrack was averaged spatially and temporally to find average pressure. The maxima and minima of spatially averaged pressure versus time were selected for each compression wave to find the average maximum pressure and average minimum pressure, respectively. The difference between average minimum and maximum pressure was calculated for each time point in the simulation and then averaged to find the difference between average minimum and maximum pressure.
Using parameter values of Womersley number Wo=0.30, flexural stiffness EI=2.3×10−12 N m2 and frequency fH=0.554 Hz, the simulation of the valveless, tubular contracting region with pericardium had a mean flow speed of 0.612 mm s−1 and a peak flow speed of 3.6 mm s−1, while the no-pericardium simulation had a mean flow speed of 0.302 mm s−1 and a peak flow speed of 4.2 mm s−1. The ratio of mean flow speed to peak flow speed was 0.17 and 0.07 for the pericardium and no pericardium cases, respectively. As with many of the simulations during parameter sweeps, the peak flow speeds of the no-pericardium simulation were higher than those of the pericardium simulation, but the mean flow speed values were twice as high for the pericardium simulation.
In addition to the higher mean flow speeds, the pericardium limited local distention of the myocardium during contraction, particularly at the junction between the myocardium and the rigid circulatory racetrack (Fig. 5). This mechanical benefit is clear in Fig. 5E, where the pericardium and myocardium touch briefly, as compared to Fig. 5B where there is nothing to prevent the myocardium from bulging out. This results in higher mean flow speeds and greater travel of the fluid markers at the same time in each simulation.
Increasing heart beat frequency increased both mean flow speeds but led to slightly lower peak flow speeds (Fig. 6A). Simulations with a pericardium did slightly better than those with no pericardium for frequencies less than 3 Hz; this difference increased at frequencies above 3 Hz. Mean-to-peak flow speed ratios were consistently better for the with-pericardium case, and the difference increased with increasing frequency (Fig. 6B).
For comparison over several orders of magnitude of Wo, only the dynamic viscosity, μ, was varied. Simulations with a pericardium outperformed simulations without a pericardium in terms of mean flow speeds for Wo (see Eqn 14 in Materials and methods) less than 1 (Fig. 6C). This difference was reduced to zero for Wo≥1. Peak flow speeds were generally lower for simulations with a pericardium than those without. Mean-to-peak speed ratios showed the largest difference around the Wo of the animal, with the pericardium simulations having ratios about 2 times higher than the no-pericardium simulations.
Within pericardium simulations, there seemed to be no relationship between myocardium flexural stiffness (EI, see Materials and methods for description) and mean flow speed. For simulations with no pericardium, mean flow speed decreased as the myocardium became more flexible while the peak flow speed remained the same (Fig. 6E). As EI increased, the ratios for the two conditions converged (Fig. 6F), with the no-pericardium case exceeding the mean flow speed of the pericardium case at values of EI>1×10−10 N m2. Mean-to-peak speed ratios for pericardium simulations were consistently higher than simulations without a pericardium.
The pericardium generally reduced the maximum average pressure and the difference between maximum and minimum pressure within the circulatory system compared with simulations without a pericardium while not affecting the average pressure (Fig. 7). When Wo was varied, simulations with no pericardium where Wo>0.2 had higher maximum pressures and lower minimum pressures, resulting in the circulatory system experiencing higher pressure differences (Fig. 7B). Pressure exhibited a similar non-linear relationship to peak flow speed with PM ratio (Fig. 7A). When a section of pericardium at the inflow end of the contracting region was left open, pressure differences were greater than those for simulations with an intact pericardium but lower than those for simulations with no pericardium (Fig. 7A).
Simulations demonstrated a non-linear relationship between PM ratio and flow speed (Fig. 8A) as well as PM ratio and pressure (Fig. 7A) when the diameter of the pericardium was altered to change the PM ratio. Both mean and peak flow speed were lower for PM ratios <2.4. Mean flow speed peaked at around a PM ratio of 2.4–2.5, whereas peak flow speed was maximized around a PM ratio of 3.0–3.5. Beyond a PM ratio of 2.5, mean flow speed dropped slightly to the PM ratio=∞ case (no pericardium).
Comparison of experiments and simulations
As peak dye speed was measured in the animal manipulations, a direct comparison can be made between peak dye speed from experiments and peak flow speed from simulations. Fig. 8B shows a plot of peak dye speed measured from experiments (Fig. 2D) overlaid with peak flow speed from simulations at various PM ratios. The curves are very similar, showing a similar increase in peak speed with increasing values of PM. The extremes of the experiments show the most obvious differences with simulations (no damage and severe damage), but these differences amount to less than 1 mm s−1 and simulation values are within the standard deviations of each experimental condition.
It is important to note a potential limitation of the comparison between experiments and simulations. Pericardium damage in the animal results in a decrease in the PM ratio due to an increase in the diameter of the myocardium. The simulations were constructed to determine how the PM ratio affects net flow for a constant diameter tube. Despite the difference in the way the PM ratio was handled, there is still reasonable agreement between experiments and simulations.
The pericardium improves blood flow driven by the myocardium
Tunicates possess a heart capable of driving blood through a resistive, semi-open circulatory system. The heart consists of a myocardium that pumps blood peristaltically and a stiff tissue layer called the pericardium which encases the myocardium. The function of the pericardium was unclear, although observations made by Anderson (1968) and Kalk (1970) indicate that pericardial damage disrupts fluid flow through the circulatory system.
Our observations made following experimental manipulation of tunicate hearts in vivo and computational simulations support the idea that the pericardium plays an important role in the myocardium's ability to contract and drive fluid flow. Experimental damage to the pericardium resulted in much lower peak dye speed (Fig. 2B,D), and computational simulations showed that adding a pericardium increases peak flow speed (Fig. 6) and lowers pressure differences (Fig. 7). The differences in fluid flow in simulations were large at the Wo and PM ratios of the animals.
There was a negative correlation between contraction ratio and PM ratio (Fig. 2C), indicating that the myocardium was unable to fully contract with a damaged or absent pericardium. Our experiments show that reduced contraction leads to lower dye speeds within the heart. The simulations show that the absence of a pericardium (where aneurysms are most pronounced) results in lower net flows compared with simulations with a pericardium.
For both the experiments and simulations, there was a range for PM ratio where blood flow was fastest. If the pericardium is too narrow compared with the myocardium diameter, both peak and mean speed are reduced, as the pericardial wall interferes with the contraction of the myocardium (Fig. 8A). If the pericardial space is too wide, there is no benefit of having the pericardium, and larger ratios approach the flow speeds of the no-pericardium simulations (PM ratio=∞). The simulations indicate that the maximum mean flow speeds occur at a PM ratio of 2.4–2.6, which corresponds to the unmanipulated PM ratios measured on animals with intact circulatory systems (Fig. 8B).
The increase in PM ratio after pericardial damage was correlated with a decrease in the contraction ratio and a drop in the blood flow produced by the heart (Fig. 2), and the computational simulations suggest that there is a physical explanation for this response. Heart beat frequency did not change with pericardium damage (Fig. 2A) or increased pressure within the heart tube. These observations circumstantially suggest heart activity is not controlled through a pressure-sensing mechanism within the myocardium alone.
The pericardium provides structural support for the myocardium
Kalk (1970) notes that narrow vessels produce a ‘considerable back pressure’ which the heart needs to overcome in order to move fluid. Yet, few structural components in the myocardium itself can resist extreme internal pressure, evidenced by the increase in myocardium diameter when the pericardium is damaged or removed. As the pericardium is very stiff, it can resist the internal pressures generated by contraction of the myocardium by limiting the total volume of the myocardial lumen and pericardial space. Results from our simulations show that the addition of a pericardium reduces differences between maximum and minimum pressure (Fig. 7). Our experimental results show that the myocardium closes the lumen less during contraction (Fig. 2C) as the diameter of the myocardium increases (reflected in the PM ratio), which in turn lowers the peak speeds produced by the heart (Fig. 2D).
It appears as if the myocardium cannot produce enough force to close the lumen sufficiently against internal pressure without the pericardium's support, nor can it resist the increase in volume of the lumen as circulatory pressure forces more blood into the myocardium. Schulze (1964) notes that each cell (with a 15:1 length to diameter ratio) has one myofibril aligned at a 45 deg angle with the length of the lumen, creating a helical spiral of myofibrils around the lumen (Kalk, 1970). When the myocardial diameter increases (as a result of damage to the pericardium), the fiber angle must also increase, which increases both the diameter of the myocardium and the length of the myofibrils (Kier, 2012). Although the force–length relationship is not known for tunicate cardiac muscle, likely the stretching of myofibrils exacerbates the inability of the myocardium to close the lumen during contraction.
In addition to promoting closure of the lumen, the pericardium also acts to limit flow extremes produced by the flexible myocardium by physically disrupting the amount of vertical displacement of the myocardium during contraction. Simulations of blood flow with no pericardium showed that the myocardium bulges to create local distention, despite the relatively low resistance posed by the racetrack circulatory system (Fig. 5). These aneurysms are particular severe at the point where the flexible heart tube meets the rigid circulatory vessel. The general effect of these aneurysms, marked by high peak flow speeds relative to the mean flow speed generated by pumping, more back flow and higher pressure differences, is that energy expended contracting the myocardium is wasted by expanding the flexible myocardium vertically, instead of driving flow inside the circulatory system.
Implications for development of the vertebrate circulatory system
Vertebrate hearts produce very high pressures to drive fluid flow through highly resistive tissues of the body. In order to produce these high pressures, the heart and vessels of the vertebrate circulatory system have an extensive network of connective tissue that resists both collapse under negative pressure and aneurysm under positive pressure. Weakening or failure of this support structure is often catastrophic (Vogel, 1994).
During evolution, the primitive cardiac tube increased in size, beat frequency and pressure-generating force to continue to move blood through an increasingly large and complex network of small vessels. By covering the myocardium in a stiff, outer layer, the mean blood speed produced by the tubular heart increased. The pericardium enhances mean blood flow under a variety of conditions that reflect increasing body size and metabolic rate. The pericardium addresses a major problem that the primitive cardiac tube faced with scaling during evolution by providing one mechanism to extend the performance space of the primitive cardiac tube without major reorganization of tissues.
Importantly, the pericardium in our simulations made the mean blood flow produced by the myocardium much less sensitive to the flexibility of the tube. This most likely benefited organisms both during ontogeny and evolution as the material properties of vertebrate cardiac tissue change with its level of organization: Young's modulus (E) varies from 5 kPa for single cardiomyocytes to 200 kPa for an isolated cardiac muscle fiber (Mathur et al., 2001). Reducing sensitivity to tissue material properties would allow animals to use a wide variety of tissue organizations to maintain blood flow during growth and evolution.
In this study, we used a 2D version of the IBM to perform numerical simulations of a pumping elastic heart tube surrounded by a stiff pericardium and compared the results with in vivo flows measured using dye visualization. Future work could incorporate measurements of the pressures within intact and damaged pericardia and heart tubes. It is also possible to obtain spatially resolved flow fields using micro-particle image velocimetry similar to what is described in Hove et al. (2003) and Vennemann et al. (2006). In terms of numerical simulations, an obvious next step is to move into three dimensions. In addition to resolving any 3D effects, this would allow the use of helically wound fibers to drive the motion of the heart tube. More detailed models of the electrophysiology and muscle mechanics could then be applied for the activation of muscle and generation of force. Finally, more accurate models of the morphology of both the heart and the circulatory system could be used to better resolve the flow fields and intra-cardiac pressure gradients. For this first study, however, we believe that the experimental and modeling techniques are sufficient to highlight the important role of the pericardium for tiny heart tubes.
MATERIALS AND METHODS
Maintenance of animals
Individuals of the tunicate species C. savignyi were shipped overnight from the Ciona Stock Center (UC Santa Barbara, CA, USA) in cooled containers and then maintained for 2–4 weeks in a recirculating, artificial seawater system at 20°C and salinity of 32–34 ppt. Animals were fed commercially available live phytoplankton mixture (Phyto Feast Live, Reef Mariculture Inc., Campbell, CA, USA) every day during the maintenance period.
Groups of two or three individuals were placed in a covered dish of seawater with menthol crystals prior to experiments until relaxed (approximately 120 min). Animals were considered completely relaxed if there was no reaction to the insertion of a sharp probe into a siphon. Individual animals were then removed from menthol-containing seawater to fresh seawater where the tunic was completely removed. Animals were then photographed with a standard centimeter ruler (which was used to measure relaxed body length using ImageJ; Rasband, 1997-2002) pinned with their atrial (excurrent) siphon to the left side of the body so the right side of the animal's body was exposed. A small opening was then cut into the visceral mass to expose the heart without damaging the major vessels leading to it. The heart was left tucked into the visceral mass and the branchial membrane was left intact to avoid pinching or otherwise constricting flow within the major vessels.
For experimental manipulations of the heart, the following conditions were used: (1) no manipulation to the pericardium (undamaged); (2) slight damage to the pericardium; and (3) severe damage to the pericardium. The pericardium was manipulated by clipping a small hole with scissors in a portion of the central pericardium at the proximal end of the myocardial loop (less than 0.5 mm2 area) for slight damage, and clipping out a large section of pericardium (over 0.5 mm2 area) for the severely damaged condition. In all treatments, the myocardium was left undamaged and still connected to the pericardium via the raphe; the intactness of the myocardium was verified by close inspection after damage was inflicted on the pericardium.
For quantitative measurements of blood flow speed through the myocardial lumen, a microliter volume syringe with a 36-gauge needle (NF36BV-2, o.d. 110 μm, World Precision Instruments, Sarasota, FL, USA) was inserted into the dorsal blood sinus. Fluorescein dye prepared in isotonic seawater (32 ppt) at 30 mg l−1 was injected into the lumen of the myocardium in a 5–7 μl bolus when the direction of myocardial contractions put the point of injection just upstream of the heart. The field was illuminated by a blue light source to increase emission of the fluorescent dye. The dye traveled through the lumen of the heart with contraction of the myocardium. The hole produced by the needle was small enough that it closed immediately after the needle was removed. This process was repeated for animals in all experimental conditions. Video was taken of myocardial contractions and a 0.5 mm ruler for scale with a Leica Fluorescence Stereo Microscope system (M205 FA, Leica Microsystems, Inc., Buffalo Grove, IL, USA) at varying frame rates with a green bandpass filter (GFP3, Leica Microsystems, Inc.).
Videos taken during experiments were used to measure specific features of myocardium and relative dye concentration. For kinematic measurements of peak dye speed in the myocardium, the position of the leading edge of the initial dye bolus was digitized in Graphclick (Arizona Software, Inc.) by hand for every other frame during the duration of the first contraction after dye injection. Two consecutive positions were used to calculate the displacement vector of the dye bolus between frames, and the time step between every other frame (based on the frame rate of the video) was used to calculate velocity. Peak speed was calculated by taking the magnitude of each velocity vector for each frame set, and the average of these speeds was calculated for all frame sets in the sample. This process was repeated three to five times for each animal, and the mean of these repetitions was used as the mean peak dye speed for each animal.
To quantify the relative contraction of the myocardium during a heart beat, the inner edge of the myocardium was digitized by hand in Graphclick during a relaxed state and during a contracted state within the period of each heart beat. The contraction ratio was calculated by dividing the contracted diameter measurement (Mc) by the relaxed diameter measurement (Mr) (Fig. 1B). This was repeated three times per animal, and the ratios were averaged for each animal. Similarly, the length of the myocardium was measured in Graphclick by selecting points on the outer edge of the myocardial tube between the two pacemaking regions of the heart.
The PM ratio was calculated by measuring in Graphclick the width of the pericardium at the proximal-most end of the myocardial loop (PD in Fig. 1B) and dividing that by the diameter of the relaxed myocardium within the myocardial loop (Mr in Fig. 1B), close to the point at which the pericardium was measured. This measurement was repeated once for each animal.
To measure the heart beat frequency, the videos were converted to grayscale image sequences (see Fig. 3 for two examples). Grayscale intensity is a measure of dye intensity, as the bandpass filter restricted light capture to wavelengths very close to the emission wavelength of fluorescein dye (520 nm) and the maximum concentration of dye within the lumen did not wash out the grayscale. Mean grayscale intensity (relative dye intensity) was measured in ImageJ over the area of the frame containing the heart loop for each video frame. Relative dye intensities were graphed in R (R Development Core Team, 2011) against time and the peaks of dye intensity for each heart beat were selected. The distance between each set of peaks was calculated, yielding the heart beat period, and the frequency was calculated as one over the period. The frequencies for each set of peaks (5–12 heart beats) were averaged to find the mean frequency for each animal.
Statistical analysis was performed with the standard statistics and car packages in R (R Development Core Team, 2011; Fox and Weisberg, 2011). Two-way ANOVA tests were used to test for significance between groups and trends and their ANOVA F-statistics and degrees of freedom are reported. Additionally, linear regression analysis was performed in R between continuous variables; these are reported with slopes, P-values associated with slope, y-intercepts (y0) and adjusted R2 values. P-values were then corrected using a Bonferroni multiple-comparison correction; corrected P-values are reported alongside F-statistics and degrees of freedom. The data were tested for outliers with a Bonferroni test for outliers in linear models.
Immersed boundary simulations
The immersed boundary method is an approach to fluid–structure interaction problems in which structures are immersed in a viscous incompressible fluid (Peskin, 2002). The immersed boundary formulation of these problems uses a Lagrangian description of the structural deformations and stresses of the immersed body and a Eulerian description of the fluid. Integral transforms with Dirac delta function kernels couple the Lagrangian and Eulerian frames. The equations are descritized; the singular delta function is replaced by a regularized version of the delta function.
The following outline describes the 2D formulation of the IBM, but the 3D extension is straightforward. For a full review of the method, please see Peskin (2002). The equations of fluid motion are given by the Navier–Stokes equations: (5) (6)
where u(x,t) is the fluid velocity, p(x,t) is the pressure, F(x,t) is the force per unit area applied to the fluid by the immersed boundary, ρ is the density of the fluid and μ is the dynamic viscosity of the fluid. The independent variables are time t and position x.
The interaction equations between the fluid and the boundary are given by: (7) (8)
where f(r,t) is the force per unit length applied by the boundary to the fluid as a function of Lagrangian position and time, δ(x) is a 2D delta function, X(r,t) gives the Cartesian coordinates at time t of the material point labeled by the Lagrangian parameter r. Eqn 8 applies force from the boundary to the fluid grid and adds an external force term, and Eqn 9 evaluates the local fluid velocity at the boundary. The boundary is then moved at the local fluid velocity, and this enforces the no-slip condition. The numerical parameters used for the simulations are given in Table 2. Unless otherwise noted, all simulations were performed on a 512×512 grid. The specific discretization used in this paper is described in Peskin and McQueen (1996).
Force equations are given for tether forces to hold the boundary in place, spring forces to resist stretching, and bending forces to resist bending. In a simple case where the boundary is to be held in place, boundary points are tethered to target points. The equation describing the force applied to the fluid by the boundary in Lagrangian coordinates is given by: (9)
where ftarg(r,t) is the force per unit length, ktarg is a stiffness coefficient and Y(r,t) is the prescribed position of the target boundary. The equations used to determine the forces due to the resistance to stretching, fstr(r,t) and fbend(r,t), are given as: (10)and (11)
kstr is a stiffness coefficient that is proportional to the resistance to stretching or compression and EI is the flexural stiffness. The total force is then the sum of these three forces: (12)
The Womersley number was calculated using the average heart beat frequency of the undamaged heart (ω=0.554 s−1), the diameter of the myocardium of the undamaged heart (Mr=0.8682 mm) and the kinematic viscosity of seawater (μ/ρ=5.8×10−5 m2 s−1) using the equation: (13)
Wo=0.3 for simulations other than Wo sweeps. For the study where we varied the Wo over several orders of magnitude, Wo was changed using only the dynamic viscosity.
Flexural stiffness for simulations other than EI sweeps (EI=2.3×10−12 N m2) was chosen based on the second moment of area (I) for a flat plate (Wainwright et al., 1982) and Young's modulus (E) for a single vertebrate cardiocyte (5 kPa) based on measurements by Mathur et al. (2001). As all boundaries in the simulation are assumed to be infinitely thin, the flexural stiffness was simply varied by changing the value of EI in Eqn 12.
We thank E. Newman-Smith and the Ciona Stock Center at University of California, Santa Barbara for animals; J. Liu, M. Ortwine, and the Zebrafish Core Facility at UNC-CH for assistance with equipment; A. Baird, H. Edwards and M. Williams for assistance with computational simulations; and two anonymous reviewers for their helpful and constructive feedback.
The authors declare no competing or financial interests.
L.D.W. and L.A.M. conceived of the question and experiments and wrote paper. L.D.W. designed, performed and analyzed data from animal experiments, and ran and analyzed data from computational simulations. L.A.M. designed and wrote computational simulations.
This research was funded by a National Science Foundation (NSF) DMS CAREER no. 1151478 (to L.A.M.) and by a NSF DMS Research and Training Grant no. 5-54990-2311 (to R. McLaughlin, R. Camassa, L.A.M., G. Forest and P. Mucha).
- © 2015. Published by The Company of Biologists Ltd