Flow within models of the vertebrate embryonic heart

ARTICLE IN PRESS Journal of Theoretical Biology 259 (2009) 449–461

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Flow within models of the vertebrate embryonic heart Arvind Santhanakrishnan, Nhi Nguyen, Jennifer G. Cox, Laura A. Miller  Department of Mathematics, The University of North Carolina, Chapel Hill, NC 27599, USA

a r t i c l e in fo

abstract

Article history: Received 2 August 2008 Received in revised form 16 March 2009 Accepted 20 April 2009 Available online 3 May 2009

Vertebrate cardiogenesis is believed to be partially regulated by fluid forces imposed by blood flow in addition to myocardial activity and other epigenetic factors. To understand the flow field within the embryonic heart, numerical simulations using the immersed boundary method were performed on a series of models that represent simplified versions of some of the early morphological stages of heart development. The results of the numerical study were validated using flow visualization experiments conducted on equivalent dynamically scaled physical models. The chamber and cardiac cushion (or valve) depths in the models were varied, and Reynolds numbers ranging from 0.01 to 1000 corresponding to the scale of the early heart tube to the adult heart were considered. The observed results showed that vortex formation within the chambers occurred for Reynolds numbers on the order of 1–10. This transition to vortical flow appears to be highly sensitive to the chamber and cushion depths within the model. These fluid dynamic events could be important to induce shear sensing at the endothelial surface layer which is thought to be a part of regulating the proper morphological development and functionality of the valves. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Heart development Heart chamber vortices Valve flow reversal Flow visualization Immersed boundary method

1. Introduction Ever since Chapman (1918) surgically removed the heart in chicken embryos reported the resultant malformation of the circulatory system nearly a hundred years ago, it has been widely recognized that the dynamics of fluid flow through living tissues is potentially connected to inducing morphological changes in the underlying vascular framework. The interplay of fluid flow in angiogenesis can be observed in pathological conditions such as during atherosclerosis, where the blood flow velocity in arteries deviates from the nominal laminar profile and flow reversal regions may be created. In the adult human heart, such atherosclerotic plaque regions restrict the arterial lumen (stenosis), therefore creating regions of high shear stress that can rupture and lead to pathological conditions including thrombosis, platelet accumulation and occlusion (Fung, 1996; Wootton and Ku, 1999). Detailed reviews of the role of fluid flow in vascular diseases may be found in articles elsewhere (Stanton, 1999; Boisseau, 2005; Patterson, 2005). The temporal characteristics of the reversed flow structure formed during the filling of the left ventricle in the diastolic phase of the human cardiac cycle have been observed to be related to proper functionality of the adult heart (Gharib et al., 2006). Most of the work in hemodynamics has been limited to the adult vasculature, however. From a developmental standpoint, fluid flows are found to influence several

 Corresponding author. Tel.: +1 919 843 8194.

E-mail address: [email protected] (L.A. Miller). 0022-5193/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2009.04.020

processes in organogenesis. Recent advances in quantitative flow visualization techniques at spatial scales on the order of several micrometers have made in vivo exploration of the fluid dynamics of the vertebrate embryonic heart possible (Hove 2004, 2006). Hove et al. (2003) experimentally showed that the frictional force imparted on vessel walls by the blood flow is important to proper morphological development of the zebrafish heart (Danio rerio). Serluca et al. (2002) showed that vascular flow in zebrafish embryos was necessary for glomerulogenesis to occur, and in those embryos with no circulation or where blood flow to the aorta was occluded, glomerular assembly was not found to occur. The flow introduced by the rotating monocilia present on the nodal surface of mice embryos has been observed to be a crucial factor in establishing the left–right asymmetry in vertebrates necessary for proper cardiovascular development (Nonaka et al., 2002; McGrath and Brueckner, 2003; Cartwright et al., 2004). Congenital heart abnormalities such as the hypoplastic left heart syndrome (HLHS) in which the left ventricle is either small or absent may primarily be triggered by improper blood flow to the developing ventricle in addition to other epigenetic factors (Gruber and Epstein, 2004). These studies suggest that research in cardiovascular development focusing on the myocardial activity and related epigenetic factors should also explore the role of fluid dynamics in the embryonic heart. It has long been noted that the endothelial cells lining the blood vessels respond to three kinds of biomechanical stimuli which include the fluid shear, fluid hydrostatic pressure, and cyclic strain (stretch) (Topper and Gimbrone, 1999). Shear stress levels as low as 0:2 dyn=cm2 can be sensed by cultured vascular

ARTICLE IN PRESS 450

A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

endothelial cells ex vivo through mechanotransduction (Olesen et al., 1988), the mechanism of which is not yet clearly understood (see the review by Weinbaum et al., 2003). Changes in the shear stress levels are sensed at the cell membrane, and these are known to trigger various genes encoding growth/transcription factors. Some examples include enodothelin-1 (ET-1) which is important during cardiovascular development and also as a vasoconstrictor, and the lung Kru¨ppel-like factor (KLF2/LKLF) which is expressed in high shear regions in the human aorta and regulates the endothelial nitric oxide synthase (NOS-3) that controls enzyme function during vasodilation (Groendendijk et al. 2005, 2007). Shear stress sensing has been proposed to occur at the extracellular component of the endothelial cells consisting of glycoproteins, also known as the glycocalyx (Reitsma et al., 2007; Jones et al., 2004). Although the structure of the endothelial glycocalyx has been the subject of numerous investigations, the flow patterns through this layer have only recently received attention (Secomb et al., 2001; Weinbaum et al., 2003; Smith et al., 2003; Leiderman et al., 2008). Experiments on cultured endothelial cells ex vivo show that compared to a baseline case with no shear conditions, the presence of flow imposed shear stress increases the rate of diffusion of macromolecules such as albumin by the glycocalyx into the endothelial cells (Yao et al., 2007). Alternately, recent in vivo investigations show that the primary cilia present in endocardial–endothelial cells function as shear stress sensors during avian embryonic heart development, especially when the levels of the stresses are relatively low (der Heiden et al., 2006). In addition, shear stress in the blood flow has been observed to influence angiogenesis via the formation of a mechanosensory complex that involves the vascular endothelial growth-factor (VEGF) receptors (Olsson et al., 2006).

1.1. Background Fig. 1 shows a simplified two-dimensional representation of the vertebrate embryonic heart at two different stages, where the events in the left half progress to those in the right half. In human embryos, the latter (rightmost) stage of heart development occurs roughly around the third week post fertilization while it is observed after one week in mouse embryos. The initial form of the

a v

Fig. 1. Flow in the vertebrate embryonic heart, drawn after Moorman et al. (2004). The left panel shows the initial form of the heart tube, with the flow being driven by peristaltic contractions dictated by the electrophysiological activity in the myocardium. Chamber formation is initiated at the next developmental stage as shown in the right panel, where a and v indicate the locations of the atria and ventricle, respectively, and alternating contractions of these regions force the blood to flow within the tube. The gray regions in the right panel denote the cardiac cushions that later develop into the valves required for maintaining unidirectional flow. At each stage of development, the fluid flow and cardiac activity are presented with increasing time, looking from top to bottom. Dashed arrows are used to denote the blood flow direction, while solid arrows show the direction of chamber contraction.

heart at the early stages can be viewed simply as a linear tube consisting of specialized myocardial cells that are capable of internal electrical activity. In the avian embryo, the electrical activity is initiated as early as the first week post fertilization in the form of a sinusoidal type ECG (Moorman et al., 2004) (corresponding to the morphology in the left panel in Fig. 1). Traditionally, the pumping of blood in the embryonic heart has been explained via peristaltic contractions that are created from the wave of depolarization in the myogenic cells. These rhythmic contractions are believed to provide the impetus necessary to move blood. The earliest electrical activity occurs in the pacemaker-like cells present in the most venous (caudal) end of the heart tube (leftmost end in Fig. 1 or closest to the atrium). The myocardial cells have varying degrees of intercellular coupling spatially, which allows differential wave conduction velocities to occur across the tube. The presence of this polarity between the caudal and cranial extremes ensures blood flows in a unidirectional manner. In general, peristalsis requires slower conduction velocities (Oð1 cm=sÞ in fish hearts) and hence poor intercellular coupling. Portions of this linear tube later locally expand into the chambers of the adult heart (ballooning), and cardiac cushion tissues are observed near the openings of the chambers (indicated by gray regions in Fig. 1). The atria develop dorsally and expand laterally, while the ventricles expand ventrally. At this ballooning stage, the ECG of the embryonic heart cycle closely resembles the adult (Moorman et al., 2004), and the chambers contract alternately to move blood across the tube. The cardiac cushions prevent backflow of blood at later development stages (as seen in the right half of Fig. 1 top to bottom) and eventually become the valves in the adult heart. The chambers are observed to have a greater degree of intercellular coupling and faster conduction velocities in comparison to the cardiac cushion regions. Detailed reviews of chamber formation in relation to genetic factors both from developmental and from evolutionary contexts can be found elsewhere (Moorman and Christoffels, 2003; Moorman et al., 2003). Interestingly, a recent study by Reckova et al. (2003) on chick embryos showed that the maturation of the conduction cells responsible for ventricular contraction depend upon the forces imparted by the blood flow. That the embryonic heart beat is not required for the purpose of nutrition but rather to aid in growth, shaping, and morphogenesis of the heart itself has been suggested by Burggren (2004) via a review of experimental work related to disruption of cardiac output in lower vertebrates. Specifically, these organisms were found to thrive entirely by diffusion alone in the embryonic stage, even in the absence of a heartbeat. Hove et al. (2003) measured the flow field within the heart tube of the zebrafish Danio rerio embryo using the erythrocytes as fluid parcel markers. By occluding the flow in both the inflow and outflow tracts, they showed that both chamber ballooning and heart looping as well as valve formation did not occur in the adult. They proposed that the blood flowing through the embryonic heart tube creates shear forces necessary for the formation and development of the heart valves. It is suggested that these flow driven forces provide a biomechanical stimulus to the endothelial surface layer, which then feed into the biochemical regulatory networks that initiate morphogenesis. While there is increasing evidence that points to the blood flow driven forces as being an important and essential factor influencing both proper cardiovascular development (embryo) and function (adult), most of the physical details of the fluid dynamics through the embryonic heart, especially at the level of shear sensing components in the endothelium, largely remain unclear. The objective of this paper is to examine the flow field within the heart tube using a combined computational and experimental approach. A series of models representing the morphological

ARTICLE IN PRESS A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

stages of the linear heart tube were employed for this study, and the effect of changing the Reynolds numbers from the embryonic case to the adult case was investigated.

2. Numerical method The immersed boundary method has been used successfully to model a variety of problems in biological fluid dynamics. Such problems usually involve the interactions between incompressible viscous fluids and deformable elastic boundaries. Some examples of biological problems that have been studied with the immersed boundary method include animal locomotion (Fauci and Peskin, 1988; Fauci, 1990; Fauci and Fogelson, 1993; Miller and Peskin 2004, 2005), cardiac blood flow (Peskin, 1972; Peskin and McQueen, 1996; McQueen and Peskin, 2000), cell deformation by shear flows (Bottino, 1998; Eggleton and Popel, 1998), the formation of blood clots by platelet aggregation (Fogelson, 1984), and wave propagation in the cochlea (Beyer, 1992; Givelberg, 1997). Some non-biological applications of the immersed boundary method include a flapping flag in two dimensions (Zhu and Peskin, 2002) and flow past a cylinder (Lai and Peskin, 2000) as a benchmark problem. For the following series of simulations, a tethered boundary version of the immersed boundary method was used to calculate the flow in the heart chambers. Basically, the channel and chamber walls are made nearly rigid by tethering each boundary point to a target point that does not move and does not interact with the fluid. The target and actual immersed boundary points are attached with virtual springs. The target stays fixed in space, and the spring attachments apply a force to the actual boundary that is proportional to the distance between corresponding points of the two boundaries. This force is then used to calculate the force applied to the fluid. Parabolic flow within the channel was prescribed upstream of the chambers by applying a force to the fluid proportional to the difference between the desired fluid motion and the actual fluid motion (see Appendix). The two-dimensional numerical simulations were constructed to be similar to the physical experiments. The outer dimensions of the computational domain were set to be 1 m  1 m. The diameter of the channel was set to 5 cm, and the maximum radius of each chamber was set to 10 cm. Although the physical distances in the simulations were not set to be equal to the experiments, the ratios of all corresponding lengths between the physical and the numerical models matched. A dynamic viscosity of

451

0:0235 kg m1 s1 (about 20 times that of water) was used in all simulations. The maximum velocity in the channel was chosen to achieve the desired Reynolds number, in order to dynamically match the flows seen in the experiments. The stiffness coefficients of the virtual springs attaching tether points to boundary points were chosen to reduce the deformation of the boundary to acceptable levels and was dependent upon Reynolds number. Since the Reynolds number was changed by varying the velocity, the forces generated are proportional to the target velocity used. Therefore, smaller stiffness coefficients provide equivalent maximum deformations at lower Reynolds numbers. Spatial and time step sizes were chosen to provide stability and convergence. The computational domain for the fluid was set to 600  600 spatial steps for all simulations. The boundary was constructed such that each boundary point was spaced 1=1200 m apart. The following section provides a detailed explanation of the model parameters and Reynolds numbers investigated.

3. Experimental arrangements 3.1. Physical models and facility description To validate the results of the numerical study, flow visualization experiments were conducted on a set of physical models designed to represent static morphological stages of the linear embryonic heart tube. Fig. 2 shows the schematic top view of the clear plexiglas models used for the experiments. The arrangement consisted of two identical halves each with dimensions 3:25 cm thick  30 cm length  8 cm height. Within each half-model, one side consisted of a curved chamber of maximum height C and three semi-elliptical cushion inserts of maximum height V, while the other side was a flat surface (see Fig. 2). The cushion inserts were made detachable to examine the effect of changing cushion depths for any single chamber dimension. The halves of the model were aligned such that the chambers faced each other along their longest dimension and were inserted in an optically accessible plexiglas flow tank with an 8 cm  8 cm cross section. Within the tank, the halves were placed so as to form a channel type arrangement through which fluid flowed internally (similar to the right half of Fig. 1), and were separated by a channel width d ¼ 1:5 cm measured along their smallest dimension. This distance d was used as the characteristic length scale toward calculating the Reynolds number Red of the fluid flow, and d was kept constant

Fig. 2. Schematic of physical model used for visualization of flow through the vertebrate embryonic heart tube. d indicates the channel width and was maintained at 1.5 cm for all the experimental results presented herein; U indicates the free stream velocity; C and V refer to the maximum depths of the chambers and cushions, respectively, the values of which are provided in Table 1.

ARTICLE IN PRESS 452

A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

throughout all of the experiments presented here. The maximum velocity was taken well upstream of the first chamber and was used as the characteristic velocity. The chamber and cushion depths were constructed as ratios of this length scale d, and the values of the different depths used are shown in Table 1. Several orders of magnitude of Reynolds numbers were examined for this study, ranging from 0.01 to 1000. As a result, the fluid medium used for the experiments had to be varied across this Red range. For the viscous end of the flows where Red 51 as well as for Red ¼ Oð100Þ, different compositions of solutions consisting of Karot brand light corn syrup and water were used, and these are shown in Table 2. For Red ¼ Oð10Þ and Oð1000Þ, water was used as the fluid medium. For Red oOð100Þ, the flow was driven by a peristaltic pump manufactured by Fisher Scientific Inc. This pump consisted of a rotary mechanism that deformed a local section of the inlet/outlet tube that emerged from the tank containing the fluid, and sent peristaltic waves at a regular interval to produce unidirectional flow. The average flow velocities imparted by this pump within the heart model upstream of the chambers were Oð1 mm=sÞ, and this changed marginally depending on the fluid used. As the peristaltic pump flow speeds were not sufficiently fast to obtain Red 4Oð10Þ, these visualization experiments were conducted in another flow tank with the same 8 cm  8 cm cross section as the flow tank described above. This tank was driven by a propeller attached to a variable speed motor (Vogel and LaBarbera, 1978) and capable of faster flow, with the average velocity on the order of 10 cm/s. To ensure uniform flow at the inflow, flow straighteners consisting of a dense array of small circular cylinders were placed against the upstream and downstream sides of the heart models in the tank.

Table 1 Table of parametric variations examined in the numerical simulations; d indicates the chamber width and was maintained at 1.5 cm in all experiments.

3.2. Diagnostic methods For the flows where Red 51 and Red  100, dye flow visualization with commercial food coloring was used to examine the flow field through the heart tube models. A small amount of the dye had to be mixed with a minor volume of the baseline fluid chosen for the particular Red prior to injecting it into the flow. This was done to balance the density of the dye solution within the time scales of the experiment. For experiments at the highly viscous end, the dye solution was used to draw straight lines upstream of the chambers in the stagnant fluid. The pump was then turned on, and these lines advected with the flow. For higher Re flows, the dye solution was injected continuously into the already moving fluid to compensate the faster dispersion of coloring. For the intermediate Red ¼ Oð1210Þ where inertia and viscosity are almost equally important, the effect of introducing the dye injection into the flow was non-negligible, and the flow field disturbances created as a result did not settle within the time of the measurement. To visualize the flow at these Reynolds numbers, the fluid (water) was mixed with a small amount of bromothymol blue (0.8 g for 4 L), which is a pH indicator that changes from yellow to a deep blue color when the pH of the solution is between 8 and 9. The solution was filtered and balanced using potassium chloride and hydrochloric acid to adjust the pH value of the fluid to just under the color transition point. A thin copper wire electrode was immersed in this solution along the central plane of the tank, upstream of the chambers to act as the cathode. Another copper wire was placed well upstream of the model near the tank inlet to act as the anode. When a DC voltage in the order of 10 V was applied across the system, electrolytic action changed the local pH of the fluid surrounding the cathode resulting in a color change of the line region to blue, which was advected with the flow when the pump was switched on. For more details on this electrode activated pH method, please see Baker (1966) and Lister (1990). A Canon EOS 5D camera with a

Table 2 Fluid media used in the flow visualization experiments.

C=d

V=d

Red

0.5 1 1 1 1.5 2 2 2

– 0.12 0.34 – – 0.12 0.34 –

0.01, 30 20 20 0.01 30 0.01 0.01 0.01, 0.1, 1, 10, 20, 30, 40, 50, 100

OðRed Þ

Composition

0.01 0.1 10 100 1000

84.375% corn syrup + 15.625% water 75% corn syrup þ25% water 100% water 6.25% corn syrup þ93:75% water 100% water

Propeller/Peristaltic Pump

flow tank

model insert

Fig. 3. Top view of the experimental setup showing the flow chamber with the physical model inserted; arrows indicate the flow direction.

ARTICLE IN PRESS A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

maximum pixel resolution of 4368 pixels  2912 pixels was used to capture an image sequence, or a 3.3 megapixel Sony Handycam was used to acquire a video of the flow field from which individual frames were extracted. Many of the experiments examined two windows of observation (one for each chamber) separately, and as a result have a certain extent of overlap between the two images. Flow visualization at Red ¼ 1 required the pH indicator method, but the streamlines were not uniform in the fluid consisting of corn syrup and water. This case will not be presented herein for the purpose of clarity.

4. Results The first part of this section examines the effect of changing Red across a constant chamber size, in the absence of cushions.

The fluid dynamic transitions observed in the various flows and their dependence on the chamber and cushion configurations will then be presented. In each case, the computations are compared to the physical experiments in order to validate the numerical model. 4.1. Effect of varying Red Figs. 4 and 5 show the flow field streamlines for the deepest chamber with no cushions from the numerical simulations and experiments for flows across six orders of magnitude of Red varying from 0.01 to 1000. This Re range corresponds to the scale of the vertebrate embryonic heart to the mature adult heart. The deepest chamber with C=d ¼ 2 was chosen as the model of interest to examine the Re changes. To isolate any other effects, no cushion inserts were used. It was expected that the deepest

Red = 0.01

Red =0.1

0.65

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0.4

0.35 0.2

0.3

0.4

0.5 Distance

0.6

0.7

0.8

0.35 0.2

0.3

0.4

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0.4

0.35 0.3

0.4

0.5 Distance

0.6

0.7

0.8

0.6

0.7

0.8

Red =30

Red =1 0.65

0.2

453

0.5 Distance

0.6

0.7

0.8

0.35

0.2

0.3

0.4

0.5 Distance

Red =100 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.2

0.3

0.4

0.5 Distance

0.6

0.7

0.8

Fig. 4. Streamlines of the flow field from the numerical simulations for Red variation on the deepest chamber ðC=d ¼ 2Þ with no cushions. Flow direction is from left to right.

ARTICLE IN PRESS 454

A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

Fig. 5. Flow visualization for varying Red on the deepest chamber ðC=d ¼ 2Þ with no cushions. Flow direction is from left to right.

chamber should initiate eddy formation first as compared to a shallow chamber. In general, the simulations agree overall to the experimental visualizations. For Red ¼ 0:01, the flow does not separate in either chamber. In other words, no closed vortices are present. The streamlines of this highly viscous flow bend around the curved chambers and do not exhibit flow reversal. In the context of classical cavity and backward-facing step flows where vortices are seen to form in the chamber at similar Re and lower aspect ratios (Dyke, 1982), the geometry is such that the flow encounters a sharp corner followed by a rectangular cavity downstream. The change in pressure gradient is made less sharp in the case of the heart chambers, and as a result the flow is able to negotiate the chambers without forming isolated recirculation regions. This is also seen in the experiments (Fig. 5(a)) where the lines of dye drawn in the flow curve around the chambers and do not form any vortices. When Red is increased to 0.1, the above trend is seen in both the computations and the experiments, and no flow reversal is observed in the chambers. At Red ¼ 1, the majority of the flow is fully attached to the model walls, but a small region appears at the bottom of both chambers with no streamlines from the main flow. The flow within this region in the chambers was observed to be substantially slower than the primary channel flow. Although not presented herein, this formation of a ‘‘dead-region’’ was also seen in the context of experiments. Increasing the Red to Oð10Þ establishes a new regime in the flow, and vortex formation is clearly observed within both the chambers in the computations in Fig. 4(d) for Red ¼ 30. This is also seen in the pH-indicator flow visualization image obtained from the experiments at Red ¼ 25 (Fig. 5(c)) where there is clearly reversed flow present in the form of a large ‘‘dead-region’’ in both the chambers. Unlike the simulations, the size of these chamber vortices is not symmetric in the experiments, with the first

chamber exhibiting nearly no forward flow compared to the second. Note that in the case of experiments where the flow field in each chamber was observed individually, the combined image showing streamlines from both chambers would not overlap perfectly, on account of both the flow and the experiment conditions. The flow within the chambers was observed from the simulations to be highly decelerated in comparison to the channel flow velocity. As a result, it was not possible to experimentally resolve the differential features of the entire flow field (in the channel and chambers). The region of reversed flow within the chambers is observed as an empty pocket with no advection of the color change as compared to the main flow. The presence of the vortex at Red ¼ Oð10Þ was verified by placing the pH indicator probe close to the chambers so as to resolve the slow moving flow region (not shown). The size of the vortices increases from being a small fraction of the chamber at Red ¼ Oð10Þ to occupying the entire chamber at Red ¼ Oð100Þ. At this Re, the chamber vortices are symmetric in both experiments and simulations. With further increase in Red , the chamber vortices become unsteady in their structure, as expected (see Fig. 5(e) for Red ¼ 1185).

4.2. Flow at low Reynolds numbers The Reynolds numbers of the flow within the embryonic heart are sufficiently low ðOð0:01  2ÞÞ owing to the high viscosity of the embryonic blood, so that the flow can be treated as fully laminar (Moorman and Christoffels, 2003). In the context of the rigid walled approximations that the numerical and physical models are herein examined, there was no occurrence of vortex formation via flow separation in the chambers. Furthermore, the fluid used

ARTICLE IN PRESS A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

in this study exhibited a linear stress–strain behavior in shearing (Newtonian), and consequently did not mimic the shear thinning, viscoelastic behavior that is characteristic of blood due to the presence of erythrocytes. To investigate if there was a possibility of inducing flow reversal via modifying the morphological aspects at this flow scale, the chamber and cushion depths were varied, and some representative results from both simulations and experiments are shown in Fig. 6. It is seen that the introduction of a shallow chamber (C=d ¼ 1, Fig. 6(a)) compared to those seen in Figs. 4(a) and 5(a) does not change the flow field, and no vortex formation or flow reversal is observed, most ostensibly due to the extremely low momentum of the flow. At these scales, inertia plays an insignificant role in the transport of the fluid, and the flow is driven by a combination of pressure gradient and viscous diffusion. As a result, the fluid tends to adhere strongly to the walls. In the event of flow separation due to a rather adverse pressure gradient, flow reattachment required for forming a closed vortex at these low Reynolds numbers would be expected to occur for lower aspect ratio chambers (length/width 51). For the deepest chamber geometry, the introduction of cushions into the flow domain also appeared to provide no means to introduce vortices in the flow field, as seen in Fig. 6(c) and (d). It was expected that placement of the cushions would present the incoming flow a steeper (more adverse) pressure gradient and thereby promote flow separation from the solid boundary. However, the flow field was not responsive to changes in flow geometry at this scale. In fact, repeating the chamber and valve morphological studies even at Red ¼ 0:1 showed no effect on the fully attached viscous flow field (not shown). Another cushion arrangement was attempted to accomplish 95% blockage (as opposed to 68% in Fig. 6(c) and (d)), and although not shown here on account of the difficulty in density matching the dye solution (to keep the dye afloat) throughout the long time scale of these experiments, there appeared to be no major impact on the flow field. Thus, within the bounds of the approximations and

Red =0.01,C/d =1

455

geometric modifications examined herein for flow conditions similar to the embryonic heart, no vortex formation was observed. To examine the flow field for other features of the linear embryonic heart tube within the scope of these models, the chambers were placed in an orientation symmetric about the centerline axis shown in Fig. 3 such that they faced each other instead of being offset (as in Fig. 5). This was done to grossly mimic the morphology of the outflow part of the embryonic heart tube (rightward to the ventricle region in Fig. 1) herein referred to as the bulbus. This anatomical structure is usually observed in the early stages of heart development in bony fishes (as bulbus arteriosus, which disappears with development) and mammals (as bulbus cordis, a portion of which later forms the right ventricle in humans), and throughout the life cycle of amphibians (as bulbus cordis). Owing to the polarity in the electrical configuration of the embryonic heart, the conduction velocities are the slowest in the region between the ventricle and the bulbus (Moorman and Christoffels, 2003). Fig. 7 shows the computational and experimental streamlines of the flow through the bulbus arrangement at Red ¼ 0:02 for the deepest chamber and cushion geometries. The flow upstream of the chambers is fully attached, even with the presence of cushions. The flow entering the chambers is attached initially, after which it separates from the solid walls and rejoins gradually at the exit of the chamber. Although not observed here due to highly decelerated flow within the chamber, it is possible that there exists two counter-rotating recirculation zones in the empty regions in the flow field.

4.3. Conditions for vortex formation In this section, the Red ¼ 10 range where chamber vortices first appear (as seen in Fig. 5) will be explored in detail. Fig. 8 shows the effect of changing Red for the deepest chamber in the absence of any cushions. The Red values were chosen around the value of

Red =0.02,C/d =1

0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.2

0.3

0.4

0.5 0.6 Distance

0.7

0.8

Red =0.01,C/d =2,V/d =0.34

Red =0.07, C/d =2, V/d =0.34

0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.2

0.3

0.4

0.5 0.6 Distance

0.7

0.8

Fig. 6. Embryonic heart type flow showing the effect of varying chamber and cushion depths: (a) and (c) present computational streamlines, (b) and (d) show the experimental flow visualization.

ARTICLE IN PRESS 456

A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

Computation.

Experiment.

0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.2

0.3

0.4

0.5 0.6 Distance

0.7

0.8

Fig. 7. Flow field streamlines for a bulbus-type arrangement: Red ¼ 0:01, C=d ¼ 2, V =d ¼ 0:34.

Red =10,computation.

Red =20,experiment.

0.65 0.6 0.55 0.5 0.45 0.4 0.35

0.2

0.3

0.4

0.5 0.6 Distance

0.7

0.8

Red =50,experiment.

Red =50,computation.

0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.2

0.3

0.4

0.5 0.6 Distance

0.7

0.8

Fig. 8. Effect of varying Red of the flow field near the vortex formation range for deepest chamber ðC=d ¼ 2Þ with no cushions.

Red  25 at which the transition was first observed in both experiments and simulations for the deepest chamber with no cushions as shown in Figs. 4 and 5. For Red ¼ 10 in the simulations, it can be seen that the transition to vortical flow is ‘‘weak’’, and the flow is very similar to what was seen in the case of Red ¼ 1 in Fig. 4(c). Increasing the value of Red to 25 results in the formation of two distinct chamber vortices (Fig. 8(c)) that nearly are of the same size as those observed for Red ¼ 100 (Fig. 4(c)). Experimental visualizations show that lowering the Red from the value of 25 required for vortex transition to 20 affects the symmetry of the process across the two chambers (Fig. 8). It can be observed in Fig. 8(b) that while the flow is fully attached in the first chamber while some reversal starts to occur in the second chamber. Increasing the Red to 50 experimentally shows the formation of a closed vortex in the second chamber and the decelerated (and hence unresolved in the visualization) reverse flow ‘‘dead region’’ in the first chamber (Fig. 8(d)). Figs. 9 and 10 examine the effect of changing the chamber dimension near the

Red required for vortex formation without cushions in these models. In general, the vortex formation process occurs first for the geometry with the deeper chamber, as expected (compare Fig. 9(a) to (b)). This is most ostensibly on account of the more adverse pressure gradient presented by the chamber with a larger depth, as discussed earlier. It is interesting to note that the vortex observed for Red 425 in Fig. 8(c) is not truly an independent condition, rather it depends on the geometry of the model used. The experimental visualizations support the trend observed in the simulations, although the test conditions are not exactly similar. The last set of experiments and simulations were conducted to determine if there was a possibility of creating flow reversal within a shallow chamber at Red close to transition in the deepest chamber, through the addition of cushions. The visualizations in Fig. 12 clearly show that by introducing the cushions at Red ¼ 35 flow, vortex formation could be induced in a rather shallow chamber geometry with C=d ¼ 1. The more the flow is blocked by the cushions, the higher the flow velocity locally through the

ARTICLE IN PRESS A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

C/d =0.5.

C/d =1.

0.65

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0.4

0.35 0.2

0.3

0.4

0.5 0.6 Distance

457

0.7

0.8

0.35 0.2

0.3

0.4

0.5 0.6 Distance

0.7

0.8

Fig. 9. Flow field streamlines from the computations showing the effect of varying C=d for Red ¼ 40 with no cushions.

Fig. 10. Flow field streamlines from the computations showing the effect of varying chamber depths at Red ¼ 25 with no cushions.

0.65

0.65 0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0.4

0.35 0.2

0.3

0.4

0.5 0.6 Distance V/d = 0.12.

0.7

0.8

0.35 0.2

0.3

0.4

0.5 0.6 Distance V/d = 0.34.

0.7

0.8

Fig. 11. Flow visualizations showing the effect of varying cushion depths for Red ¼ 20 on the deepest chamber C=d ¼ 2.

Fig. 12. Flow visualizations showing the effect of varying cushion depths for Red ¼ 35 and chamber depth C=d ¼ 1.

contracted region where it is subjected to an adverse pressure gradient due to the curvature of the cushions, and the larger the chamber vortices as seen from comparing Fig. 12(a) and (b). The simulations shown in Fig. 11 establish the same result seen in the experiments.

5. Summary The conditions required for vortex formation are significantly affected by flow Reynolds number and are highly sensitive to the chamber and cushion dimensions within the band of Re where

ARTICLE IN PRESS 458

A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

this occurs. In the context of heart development, this implies that the morphology of the developing heart is important to the dynamics of fluid flow within. As the early chambers bulge from the heart tube, flow velocities may be reduced near the chamber walls. If chamber vortices form, flow reversals can occur near the cardiac walls, changing the magnitude and direction of shear stress. As the cardiac cushions begin to form, flow velocities and shear stresses in the atrio-ventricular canal can increase dramatically. If the formation of the cushions induces vortex formation within the cardiac chambers, then the magnitude and direction of the wall shear stress can change significantly. The dependence of the flow upon Reynolds number, chamber depth, and cushion height renders the problem not universally scalable (say even within different stages in the embryonic heart maturation). As a result, the blood flow at each developmental stage of the heart should be carefully investigated. 5.1. Comparison to in vivo and in vitro studies For Reynolds numbers on the order of the embryonic heart tube, vortices were not observed in either the flow through the physical models or the numerical simulations. The in vitro flow visualization of the zebrafish embryonic heart by Hove et al. (2003), however, revealed the formation of a fluidic jet-like flow field when the blood is pumped from the ventricle to the bulbus, with a counter-rotating vortex ring being formed not on account of flow separation but rather due to the complete blockage of the bulbus exit at this stage of the cardiac cycle. These measurements relied on employing the red blood cells as passive fluid markers. Interestingly, on the scale of the embryonic heart, the diameters of these cells are on the same order of magnitude as the cushion entrances. Also, the red blood cells are flexible and unevenly shaped. As a result, the vortices observed in their study might be affected by the tumbling and rotation of the red blood cells, although the effects of this might be rather non-uniform in time. Vennemann et al. (2006) performed in vitro flow field measurements on a chicken embryo using liposome particles as tracers that were artificially introduced into the flow, and these particle had length scales that were an order of magnitude smaller compared to the red blood cells. They synchronized their measurements to different portions of the cardiac cycle, and did not report any vortex formation events. The importance of vortex formation stems from the shear sensing capability of the endothelial cells which can then feed into the biochemical pathways to activate genes required for morphological development. In the presence of a vortex region, flow reversal is observed, and it is not only the magnitude but the gradient of shear stresses and direction of flow that are significant factors to the shear sensing by the endothelial cells (Jones, 2005). The static approximation of the actual pumping dynamics in the current study, coupled with the lack of particulate content (as in red blood cells in the actual embryonic heart) and absence of a pulsatile input flow, is thought to all play an important role toward why the jet-like flow field was not observed in the physical model experiments. 5.2. Comparison to other numerical studies of the embryonic heart A couple of recent numerical investigations have described the blood flow through sophisticated three-dimensional models of the vertebrate embryonic heart. DeGroff et al. (2003) used a sequence of two-dimensional cross-sectional images to reconstruct the three-dimensional surface of human heart embryos at stages 10 and 11. In their paper, the heart walls did not move, and steady and pulsatile flows were obtained using finite volume CFD. Similar

to our results, they did not report the presence of chamber vortices. Their study also showed that streaming was present in the heart tube (particles released on one side of the lumen did not cross over or mix with particles released from the opposite side). Liu et al., 2007 quantified the hemodynamic forces on a three-dimensional model of a chick embryonic heart using a finite element model. They focused on pulsatile flow through the outflow tract during stage HH21 (after about 3.5 days of incubation) and included flexibility in the walls of the tract. They did not include cardiac cushions in their simulations. Similar to our results, maximum velocities were observed in regions of constrictions. Vortices were observed during the ejection phase near the inner curvature of the outflow tract. The maximum Reynolds number used in their simulations was 6.9, which is close to point where we observed the transition to vortical flow in our simplified models. It would be interesting to explore the effect of pulsatile flow and wall flexibility on the formation of vortices in the embryonic heart. 5.3. Limitations of the model In order to explore a large parameter space, a number of simplifications to the physical and numerical models were made that limit the degree to which the results can be interpreted in the context of the actual embryonic heart. Our goal in this paper was to report fluid dynamic transitions in simplified models due to changes in morphology and Reynolds number that can then be rigorously explored with three-dimensional simulations and sophisticated physical models. Two of the limitations of the model are that the flow is steady and that the heart walls are rigid. It is certainly possible that unsteady effects and heart wall flexibility may change the point at which chamber vortices form which would also change the direction and magnitude of shear felt at the heart walls. Another limitation of the model is that the flow is two dimensional, and blood flow in the embryonic heart is inherently three dimensional. The three dimensionality of the flow is clearly observed in the CFD results of DeGroff et al. (2003) and Liu et al. (2007) as fluid particles wind through the complicated morphology of the heart tube. Finally, we expect that the actual contraction and expansion of the cardiac chambers will alter the intracardiac fluid dynamics and shear stresses felt along the heart wall. Extending our results to three dimensions and incorporating the actual contraction and expansion of the heart chambers will be the focus of future work. 5.4. Applications to other biological problems The results of the study may also be important to several anatomical structures that are found in humans that closely resemble the geometry of the physical model employed herein. The aortic sinus found in the ascending aorta, for example, has a set of three semilunar valves all of which resemble the chambers in the physical models. It is well known that the vortex structure formed in the aortic sinus is necessary for proper functionality of the aortic valve closure and to prevent backward flow of blood into the ventricle. Peskin and Wolfe (1978) reviewed this phenomenon briefly, and presented an analytical and numerical study of valve closure. Peacock (1990) examined the structure of the vortex experimentally and observed that localized turbulence might be introduced depending on the physiological condition. The chamber vortices in the developing heart could potentially be important in assisting the cardiac cushion closure (as depicted in Fig. 1). It is possible that the flow within the sinus is important to proper morphological development of the aortic valves during later stages of cardiovascular system development in humans.

ARTICLE IN PRESS A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

Another possible application of flows at this Re relates to flow across different portions of the adult human pulmonary system (McDonald, 1997). These flows encompass a wide range of Reynolds numbers spanning from roughly 0.01 in the ciliary regions to 1000 and greater in the trachea. The alveolar sacs and bronchi anatomically show the presence of surface roughness and pocket-like regions, and the physical models constructed herein could be used to understand the fluid flow within these morphological features.

Acknowledgments The authors would like to thank the Applied Math Fluids Group at UNC-CH for their valuable insight into the implications of this work and its relation to classical problems. This work was supported by a Burroughs Wellcome Fund CASI award ID # 1005782.01.

459

The force equations are given by f ðr; tÞ ¼ ktarg ðYðr; tÞ  Xðr; tÞ

(5)

F ext ðx; tÞ ¼ kext ðuðx; tÞ  utarg ðx; tÞ 2   !3 0:5  x 2 7 6 U max 1  7 d=2 utarg ðx; tÞ ¼ 6 5 4 0

(6) (7)

Eq. (5) describes the force applied to the fluid by the boundary in Lagrangian coordinates. f ðr; tÞ is the force per unit area, ktarg is a stiffness coefficient, and Yðr; tÞ is the prescribed position of the target boundary. Eq. (6) describes the force applied to the fluid in Eulerian coordinates as a result of the external force necessary to drive the channel flow. Eq. (7) describes utarg , the desired parabolic flow profile upstream of the chambers where d is the channel diameter and U max is the desired maximum velocity in the channel. Discretization

Appendix The immersed boundary method provides a way to solve coupled fluid-structure problems based on forces. The version of the immersed boundary method used in this paper has been described in detail by Peskin and McQueen (1996). A review of the numerical method follows below, and more details on the immersed boundary formulation can be found in Peskin (2002). The equations of motion for the fluid are given by the Navier–Stokes equations in Eulerian form:   @uðx; tÞ r (1) þ uðx; tÞ  ruðx; tÞ ¼ rpðx; tÞ þ mDuðx; tÞ þ Fðx; tÞ @t

r  uðx; tÞ ¼ 0

(2)

where uðx; tÞ is the fluid velocity, pðx; tÞ is the pressure, Fðx; tÞ is the force per unit volume applied to the fluid by the immersed boundary, r is the density of the fluid, and m is the dynamic viscosity of the fluid. The independent variables are the time t and the position x. Eq. (2) is the condition that the fluid is incompressible. We use three-dimensional terminology in this description even though our model is two dimensional. This is consistent if we regard the model as a cross section of a threedimensional apparatus that looks the same in every cross section parallel to the x, y plane. The interaction equations between the fluid and the boundary are given by Z (3) Fðx; tÞ ¼ f ðr; tÞdðx  Xðr; tÞÞ dr þ F ext ðx; tÞ Z @Xðr; tÞ (4) ¼ UðXðr; tÞ ¼ uðx; tÞdðx  Xðr; tÞÞ dx @t where f ðr; tÞ is the force per unit area applied by the boundary to the fluid as a function of Lagrangian position and time, F ext is the external force applied to drive the flow, dðxÞ is a two-dimensional delta function, Xðr; tÞ gives the Cartesian coordinates at time t of the material point labeled by the Lagrangian parameter r. Eq. (3) applies force from the boundary to the fluid grid and adds an external force term, and Eq. (4) 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. Each of these equations involves a two-dimensional Dirac delta function d, which acts in each case as the kernel of an integral transformation. These equations convert Lagrangian variables to Eulerian variables and vice versa.

The system of differentio-integral equations given by Eqs. (3)–(7) was solved on a rectangular grid with periodic boundary conditions in both directions. The Navier–Stokes equations were discretized on a fixed Eulerian grid, and the immersed boundaries were discretized on a Lagrangian array of points. To begin, consider the discretization of the equations that describe the force applied to the fluid as a result of the distance between the actual boundary and the tether points (Eq. (5)). Let Dt be the duration of the time step, let Ds be the length of a boundary spatial step, let n be the time step index, and let s define the location of a boundary point in the Lagrangian framework (s ¼ mDs, m is an integer). Also, let X n ðsÞ ¼ Xðs; nDtÞ, Y n ðsÞ ¼ Yðs; nDtÞ, un ðxÞ ¼ uðx; nDtÞ, and pn ðxÞ ¼ pðx; nDtÞ. Using these definitions, the force applied to the fluid by the boundary, f , defined by Eq. (5), is discretized as follows: n

f ¼ ktarg ðY n ðsÞ  X n ðsÞÞ

(8)

The external force applied to the fluid in order to drive the flow, F ext , defined in Eq. (6) and the desired flow profile upstream, utarg ðx; tÞ, defined by Eq. (7) can be discretized as follows: ( kext ðun  untarg Þ; 0:5  d=2  jh  0:5 þ d=2; 10  i  15 n (9) F ¼ 0 otherwise 2  2 ! 3 0:5  jh 7 6 U max 1  7 d=2 (10) untarg ði; jÞ ¼ 6 5 4 0 where x ¼ ðih; jhÞ, i and j are the integers, and h is the mesh width. Since the Lagrangian array of boundary points does not necessarily coincide with the fixed Eulerian lattice used for the computation of the fluid velocities, a smoothed approximation to the Dirac delta function (Eqs. (3) and (4)) is used to handle the fluid–boundary interaction. This approximate delta function applies a force from the boundary to the fluid and interpolates the fluid velocity. The approximate two-dimensional Dirac delta function, dh , used in these calculations is given by the following equations:  x  y dh ðxÞ ¼ h2 f f (11) h h 8     < 1 1 þ cos pr if jrj  2 2 (12) fðrÞ ¼ 4 : 0 otherwise where h is the size of the spatial step in the Eulerian grid and x ¼ ðx; yÞ. Peskin and McQueen (1996) have previously described

ARTICLE IN PRESS 460

A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

this choice of the delta function. Other approximate delta functions have also been used successfully (see Peskin, 2002). Using this choice of dh , Eqs. (3) and (4) can be discretized as follows: X n F n ðxÞ ¼ f ðsÞdh ðx  X n ðsÞDs þ F next ðxÞ (13) s

U nþ1 ðsÞ ¼

X

2

unþ1 ðxÞdh ðx  X n ðsÞh

(14)

x n

X nþ1 ðsÞ  X ðsÞ ¼ U nþ1 ðsÞ (15) Dt P where s is the sum over all of the discrete collection of points of P the form s ¼ mDs, and x is the sum over of the discrete points of the form x ¼ ðih; jhÞ, where i, j and m are the integers and h is the mesh width. The final step is to discretize the incompressible Navier–Stokes equations. There are many different schemes that have been used here as part of the immersed boundary method (e.g., Peskin and McQueen, 1996; Lai and Peskin, 2000; Peskin, 2002). The method used here is implicitly defined as follows:  nþ1  2 X u  un  nþ1 þ Sh ðun Þun þ D0 pnþ1 ¼ m Dþ þ Fn a Da u Dt a¼1

r

(16)

D0  unþ1

(17)

where r is the density of the fluid and m is the viscosity of the fluid. D0 is the central difference approximation to r. It is defined as follows: D0 ¼ ðD01 ; D02 Þ ðD0a fÞðxÞ ¼

(18)

fðx þ hea Þ  fðx  hea Þ

(19)

2h

where fe1; e2g is the standard basis of R2 . The operators D a are forward and backward difference approximations of @=@xa . They are defined as follows: ðDþ a fÞðxÞ ¼

fðx þ hea Þ  fðxÞ

h fðxÞ  fðx  hea Þ  ðDa fÞðxÞ ¼ h

(20) (21) P2

þ



Thus, the viscous term given as a¼1 Da Da is a difference approximation to the Laplace operator. Finally, Sh ðuÞ is a skewsymmetric difference operator which serves as a difference approximation of the nonlinear term u  ru. This skew-symmetric difference operator is defined as follows: Sh ðuÞf ¼ 12u  D0h f þ 12D0h  ðufÞ

(22)

For further discussion, see Peskin (2002). Since the equations are linear in the unknowns unþ1 and pnþ1 , the fast Fourier transform (FFT) algorithm was used to solve for unþ1 and pnþ1 from un , pn , and F n (see Press et al., 1992; Peskin and McQueen, 1996). References Baker, D.J., 1966. A technique for the precise measurement of small fluid viscosities. J. Fluid Mech. 26, 573–575. Beyer, R.P., 1992. A computational model of the cochlea using the immersed boundary method. J. Comput. Phys. 98, 145–162. Boisseau, M.R., 2005. Roles of mechanical blood forces in vascular diseases. Clin. Overview Clin. Hemorheol. Microcirc. 33, 201–207. Bottino, D.C., 1998. Modeling viscoelastic networks and cell deformation in the context of the immersed boundary method. J. Comput. Phys. 147, 86–113. Burggren, W.W., 2004. What is the purpose of the embryonic heart beat? Or how facts can ultimately prevail over physiological dogma. Physiol. Biochem. Zool. 77, 333–345. Cartwright, J.H.E., Piro, O., Tuval, I., 2004. Fluid-dynamical basis of the embryonic development of left–right asymmetry in vertebrates. PNAS 101, 7234–7239. Chapman, W.B., 1918. The effect of the heart-beat upon the development of the vascular system in the chick. Am. J. Ant. 23, 175–203.

DeGroff, C.G., Thornburg, B.L., Pentecost, J.O., Thornburg, K.L., Gharib, M., Sahn, D.J., Baptista, A., 2003. Flow in the early embryonic human heart: a numerical study. Pediatr. Cardiol. 24, 375–380. der Heiden, K.V., Groenendijk, B.C.W., Hierck, B.P., Hogers, B., Koerten, H.K., Mommaas, A.M., de Groot, A.C.G., Poelmann, R.E., 2006. Monocilia on chicken embryonic endocardium in low shear stress areas. Dev. Dynam. 235, 19–28. Dyke, M.V., 1982. An Album of Fluid Motion. The Parabolic Press, Stanford, CA. Eggleton, C.D., Popel, A.S., 1998. Large deformations of red blood cell ghosts in a simple shear flow. Phys. Fluids 10, 1834–1845. Fauci, L.J., 1990. Interaction of oscillating filaments—a computational study. J. Comput. Phys. 86, 294–313. Fauci, L.J., Fogelson, A.L., 1993. Truncated Newton methods and the modeling of complex immersed elastic structures. Comm. Pure Appl. Math. 46, 787–818. Fauci, L.J., Peskin, C.S., 1988. A computational model of aquatic animal locomotion. J. Comput. Phys. 77, 85–108. Fogelson, A.L., 1984. A mathematical model and numerical method for studying for studying platelet adhesion and aggregation during blood clotting. J. Comput. Phys. 56, 111–134. Fung, Y.C., 1996. Biomechanics: Circulation, second ed. Springer, New York. Gharib, M., Rambod, E., Kheradvar, A., Sahn, D.J., Dabiri, J.O., 2006. Optimal vortex formation as an index of cardiac health. PNAS 103, 6305–6308. Givelberg, E., 1997. Modeling elastic shells immersed in fluid. Ph.D. Thesis, Courant Institute of Mathematical Sciences, New York University, New York, New York. Groenendijk, B.C.W., der Heiden, K.V., Hierck, B.P., Poelmann, R.E., 2007. The role of shear stress on ET-1, KLF2, and NOS-3 expression in the developing cardiovascular system of chicken embryos in a venous ligation model. Physiology 22, 380–389. Groenendijk, B.C.W., Hierck, B.P., Vrolijk, J., Baiker, M., Pourquie, M.J.B.M., de Groot, A.C.G., Poelmann, R.E., 2005. Changes in shear stress-related gene expression after experimentally altered venous return in the chicken embryo. Circ. Res. 96, 1291–1298. Gruber, P.J., Epstein, J.A., 2004. Development gone awry—congenital heart disease. Circ. Res. 94, 273–283. Hove, J.R., 2004. In vivo biofluid imaging in the developing Zebrafish. Birth Defects Res. 72, 277–289. Hove, J.R., 2006. Quantifying cardiovascular flow dynamics during early development. Pediat. Res. 60, 6–13. Hove, J.R., Ko¨ster, R.W., Forouhar, A.S., Bolton, G.A., Fraser, S.E., Gharib, M., 2003. Intracardiac fluid forces are an essential epigenetic factor for embryonic cardiogenesis. Nature 421, 172–177. Jones, E.A.V., 2005. Blood flow and the mammalian embryo. Ph.D. Thesis, California Institute of Technology, Pasadena, CA. Jones, E.A.V., Baron, M.H., Fraser, S.E., Dickinson, M.E., 2004. Measuring hemodynamics during development. Am. J. Physiol. Heart Circ. Physiol. 287, H1561–H1569. Lai, M.-C., Peskin, C.S., 2000. An immersed boundary method with formal second order accuracy and reduced numerical viscosity. J. Comput. Phys. 160, 705–719. Leiderman, K.M., Miller, L.A., Fogelson, A.L., 2008. The effects of spatial inhomogeneities on flow through the endothelial surface layer. J. Theor. Biol. 252(2), 313–325. Lister, C.R.B., 1990. An explanation for the multivalued heat transport found experimentally for convection in a porous medium. J. Fluid Mech. 214, 287–320. Liu, A., Rugonyi, S., Pentecost, J.O., Thornburg, K.L., 2007. Finite element modeling of blood flow-induced mechanical forces in the outflow tract of chick embryonic hearts. Comput. Struct. 85, 727–738. McDonald, J.A., 1997. Lung Growth and Development (Lung Biology in Health and Disease), Informa Healthcare, New York. McGrath, J., Brueckner, M., 2003. Cilia are at the heart of vertebrate left–right asymmetry. Curr. Opin. Genet. Dev. 13, 385–392. McQueen, D.M., Peskin, C.S., 2000. A three-dimensional computer model of the human heart for studying cardiac fluid dynamics. Comput. Graph. 34, 56–60. Miller, L.A., Peskin, C.S., 2004. When vortices stick: an aerodynamic transition in tiny insect flight. J. Exp. Biol. 207, 3073–3088. Miller, L.A., Peskin, C.S., 2005. A computational fluid dynamics study of ‘clap and fling’ in the smallest insects. J. Exp. Biol. 208, 195–212. Moorman, A., Webb, S., Brown, N.A., Lamers, W., Anderson, R.H., 2003. Development of the heart: (1) formation of the cardiac chambers and arterial trunks. Heart 89, 806–814. Moorman, A.F.M., Christoffels, V.M., 2003. Cardiac chamber formation: development, genes, and evolution. Physiol. Rev. 83, 1223–1267. Moorman, A.F.M., Soufan, A.T., Hagoort, J., Boer, P.A.J.D., Christoffels, V.M., 2004. Development of the building plan of the heart. Ann. NY Acad. Sci. 1015, 171–181. Nonaka, S., Shiratori, H., Saijoh, Y., Hamada, H., 2002. Determination of left–right patterning of the mouse embryo by artificial nodal flow. Nature 418, 96–99. Olesen, S.P., Clapham, D.E., Davies, P.F., 1988. Hemodynamic shear-stress activates a Kþ current in vascular endothelial cells. Nature 331, 168–170. Olsson, A.-K., Dimberg, A., Kreuger, J., Claesson-Welsh, L., 2006. VEGF receptor signalling—in control of vascular function. Mol. Cell Biol. 7, 359–371. Patterson, C., 2005. Even flow: shear cues vascular development. Arterioscler. Thromb. Vasc. Biol. 25, 1761–1762. Peacock, J.A., 1990. An in vitro study of the onset of turbulence in the sinus of valsalva. Circ. Res. 67, 448–460. Peskin, C.S., 1972. Flow patterns around heart valves: a numerical method. J. Comput. Phys. 10, 252–271.

ARTICLE IN PRESS A. Santhanakrishnan et al. / Journal of Theoretical Biology 259 (2009) 449–461

Peskin, C.S., 2002. The immersed boundary method. Acta Numerica. 11, 479–517. Peskin, C.S., McQueen, D.M., 1996. Fluid dynamics of the heart and its valves. In: Othmer, H.G., Adler, F.R., Lewis, M.A., Dallon, J.C. (Eds.), Case Studies in Mathematical Modeling: Ecology, Physiology, and Cell Biology, second ed. Prentice-Hall, New Jersey. Peskin, C.S., Wolfe, A.W., 1978. The aortic sinus vortex. Fed. Proc. 37, 2784–2792. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in Fortran 77: The Art of Scientific Computing, second ed. Cambridge University Press, UK. Reckova, M., Rosengarten, C., deAlmeida, A., Stanley, C.P., Wessels, A., Gourdie, R.G., Thompson, R.P., Sedmera, D., 2003. Hemodynamics is a key epigenetic factor in development of the cardiac conduction system. Circ. Res. 93, 77–85. Reitsma, S., Slaaf, D.W., Vink, H., van Zandvoort, M.A.M.J., oude Egbrink, M.G.A., 2007. The endothelial glycocalyx: composition, functions, and visualization. Eur. J. Physiol. 454, 345–359. Secomb, T., Hsu, R., Pries, A., 2001. Effect of the endothelial surface layer on transmission of fluid shear stress to endothelial cells. J. Biorheol. 38, 143–150. Serluca, F.C., Drummond, I.A., Fishman, M.C., 2002. Endothelial signaling in kidney morphogenesis: a role for hemodynamic forces. Curr. Biol. 12, 492–497. Smith, M., Long, D., Damiano, E., Ley, K., 2003. Near-wall micro-PIV reveals a hydrodynamically relevant endothelial surface layer in venules in vivo. Biophys. J. 85, 637–645.

461

Stanton, A.V., 1999. Haemodynamics, wall mechanics and atheroma: a clinician’s perspective, Proceedings of the Institution of Mechanical Engineers. Part H: J. Eng. Med. 213, 385–390. Topper, J.N., Gimbrone, M.A., 1999. Blood flow and vascular gene expression: fluid shear stress as a modulator of endothelial phenotype. Mol. Med. Today 5, 40–46. Vennemann, P., Kiger, K.T., Lindken, R., Groenendijk, B.C.W., de Vos, S.S., ten Hagen, T.L.M., Ursem, N.T.C., Poelmann, R.E., Westerweel, J., Hierck, B.P., 2006. In vivo micro particle image velocimetry measurements of blood–plasma in the embryonic avian heart. J. Biomech. 39, 1191–1200. Vogel, S., LaBarbera, M., 1978. Simple flow tanks for research and training. Bioscience 28, 638–643. Weinbaum, S., Zhang, X., Han, Y., Vink, H., Cowin, S.C., 2003. Mechanotransduction and flow across the endothelial glycocalyx. PNAS 100, 7988–7995. Wootton, D.M., Ku, D.N., 1999. Fluid mechanics of vascular systems, diseases, and thrombosis. Annu. Rev. Biomed. Eng. 1, 299–329. Yao, Y., Rabodzey, A., Dewey, C.F., 2007. Glycocalyx modulates the motility and proliferative response of vascular endothelium to fluid shear stress. Am. J. Physiol. Heart Circ. Physiol. 293, H1023–H1030. Zhu, L., Peskin, C.S., 2002. Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method. J. Comput. Phys. 179, 452–468.