Flow dynamics control the location of sprouting and direct elongation during developmental angiogenesis

The vasculature is the first organ to form during development and it adapts to the onset of blood flow in part by forming new blood vessels through angiogenesis. The ability to control angiogenesis has important therapeutic applications for many diseases. Anti-angiogenic treatments used to treat tumours lead to transient vessel normalisation, which causes a reduction in tumour invasiveness and better delivery of chemotherapeutic drugs (Carmeliet and Jain, 2011). In wet macular degeneration, defective angiogenesis leads to the formation of leaky vessels that damage the eye and lead to vision loss. Anti-angiogenic therapies can effectively stop and sometimes reverse the progression of macular degeneration (Tah et al., 2015). As such, there is significant interest in understanding the signals that induce or inhibit angiogenesis.

In the adult vasculature, endothelial cells exist in a quiescent state wherein little proliferation occurs (Hanahan and Folkman, 1996; Hobson and Denekamp, 1984). Normal physiological levels of shear stress maintain endothelial cells in a quiescent state by inducing cell cycle arrest (Lin et al., 2000) and strengthening tight junctions (Lin et al., 2000). For sprouting to occur, however, endothelial cells must re-enter the cell cycle (Bai et al., 2014), loosen their attachment to neighbouring cells (Lampugnani and Dejana, 2007), degrade the local basement membrane (Siefert and Sarkar, 2012) and then invade the tissue. Although shear stress is known to affect many of these pathways, very little is known about the effects of shear stress on angiogenesis. Two groups recently developed micro-machined flow cells that enabled them to alter shear stress and interstitial flow as endothelial cells sprout (Galie et al., 2014; Song and Munn, 2011). Both groups found that interstitial flow enhanced angiogenesis, however, one found that shear stress inhibited sprouting (Song and Munn, 2011), whereas the other found that it induced sprouting (Galie et al., 2014). Using an embryonic model, our group has demonstrated that increasing blood viscosity, which increases shear stress, results in a decrease in vascular density (Chouinard-Pelletier et al., 2013). We showed by time-lapse microscopy that this occurs because of a decrease in sprouting angiogenesis. Changing viscosity, however, also altered other hemodynamic parameters such as the pressure drop through the network. These studies establish a role for flow dynamics in angiogenesis but a better understanding of the hemodynamic environment present in vivo during sprouting is needed to identify the causal stimuli.

The vascular network must not only deliver oxygen and nutrients to tissues, but it must also be hemodynamically efficient. If a new vessel forms in a hemodynamically inefficient location, then a high resistance to flow will be present and limited flow will pass through that vessel. The production of a hemodynamically efficient vascular network can occur either by excessive production of blood vessels followed by the pruning of inefficient vessels, or by sprouting of new vessels in hemodynamically efficient locations in the first place. Though vascular networks do overproduce vessels and prune inefficient ones (Pries and Secomb, 2014), this does not negate the possibility of a bias at the site of sprout initiation in perfused vascular networks. It therefore makes physiological sense that the forces created by blood flow could influence the angiogenic process.

To study the role that flow dynamics play in the process of angiogenesis, we used a technique to simultaneously image vessel morphology and flow dynamics by time-lapse microscopy in the capillary plexus of avian embryos undergoing vascular remodelling (Ghaffari, et al., 2015). We injected Alexa Fluor 488-labelled acetylated low-density lipoprotein (AF488-AcLDL), which specifically labels endothelial cells and macrophages (Brown et al., 1980; Hallmann et al., 1987). AF488-AcLDL labels not only perfused vessels, but also sprouting endothelial cells. We then injected red fluorescent microspheres to follow blood flow dynamics and imaged the microspheres with a high-speed camera. We analysed the pattern of blood flow present during sprouting in the arterial plexus of three quail embryos as well as the venous plexus of three additional quail embryos, for a total of 6 time-lapse recordings (Fig. 1A). Early in vascular development, a spatial separation exists such that the rostral capillary plexus is venous and the caudal plexus is arterial. The motion of the microspheres was analysed using a technique called micro-particle image velocimetry (µPIV) to calculate the velocity of the blood. Micro-PIV was only used to make velocity measurements in straight vessel segments that led in and out of a region of interest. These measurements were then used as the input for a computational fluid dynamic solver to calculate flow in the entire region of interest, around the growing sprout. Flow dynamics were analysed starting 30 min before the sprout was visible and continued until the new vessel lumenised and carried flow (Fig. 1B, yellow arrowheads). The computational solver calculates blood velocity in the entire region of interest from the μPIV inputs (Fig. 1C,D), as well as calculating shear stress, pressure drop and vorticity (Movies 1-3).

We first investigated the relationship between flow parameters and the location of sprout initiation. We found that sprouts consistently formed from vessels at lower pressure towards vessels at higher pressure (Fig. 2, n=6/6 embryos, P=0.014). The contour plots of relative pressure 30 min before sprout formation are presented, both for the arterial (Fig. 2A) and the venous (Fig. 2C) regions of the plexus. The location and direction of sprouting are identified by the black arrows (Fig. 2A,C) as observed in the green images of vascular morphology after the sprout has formed (Fig. 2B,D, yellow arrowheads). The average pressure difference between the initiation site and the opposite vessel (i.e. the one towards which the sprout is extended) was of the same order of magnitude in all embryos. Pressure differences were normalised by the distance between the two vessels. The average value was +2.1 Pa/mm (s.e.m.=0.09, n=3) for arteries and +1.7 Pa/mm (s.e.m.=0.09, n=3) for veins. Though sprouts formed from lower towards higher pressure vessels, the location and direction of sprouting was not at the site of the maximum pressure difference in any of the time-lapse experiments.

Shear stress has been shown to play a significant role in controlling angiogenesis. We therefore investigated the level of shear stress on the endothelium when sprouting initiated (i.e. 30 min prior to sprout appearance, Fig. 3). Shear stress contour plots just prior to sprouting for the three arterial embryos are shown in Fig. 3A, whereas Fig. 3B shows the level of shear stress with respect to distance along the endothelium. The black lines and arrows in Fig. 3A indicate the arc length represented in the x-axis of Fig. 3B (x=0 at the bottom of the arrow). Shear stress levels were plotted for a subsection of the lower pressure region; however plots of the entire region are available in Fig. S1. We found that for all arterial time-lapse recordings, the location of sprouting could be identified by a local minimum in shear stress (Fig. 3B, sprouting location indicated by black arrows, n=3/3 embryos). Our previous results had indicated that shear stress inhibited sprouting angiogenesis in the venous but not the arterial plexus (Chouinard-Pelletier et al., 2013). In our venous time-lapse experiments, the sprout location was not associated with any specific pattern or level of shear stress (data not shown). One occurred at a minimum; the second at a midpoint; and the third was very close to the point with highest level of shear stress in the region. Therefore, in veins, our current results do not support a correlation between the level of shear stress and the sprouting location.

Though sprouts in the arterial plexus formed at a local minimum in shear stress, this was not always the absolute minimum. We investigated flow patterns at all shear stress minima (Fig. 4A, blue arrow represents sprouting location, orange arrow represents the second shear stress minimum in Embryo 1). We observed that non-sprouting shear stress minima occurred in regions where two streams of blood flow merged (Fig. 4A, n=3/3 embryos). Vorticity is the angular velocity of a fluid, and convergence or divergence of streams results in a sign change in the vorticity. Regions where streams converge can clearly be identified in the contour plots of the vorticity because the vorticity changes from a positive to a negative value (Fig. 4B, orange arrow indicating region with a change from red to blue). Though there is significant vorticity in this region, this occurs because of the angular velocity of the fluid that results from the streams converging and not because of the presence of recirculation. We plotted shear stress (Fig. 4C,E,G) and vorticity (Fig. 4D,F,H) for the entire region at lower pressure. If points where the vorticity passes through zero are excluded from the analysis, then the lowest minimum in shear stress defines the location of sprouting (n=3/3 embryos, P=0.010).

We next tested whether we could use these observations to predict sprouting location. We analysed flow dynamics prior to sprouting in the arterial region of four new embryos (Fig. 5). The region and time point that were analysed were chosen by one author but analysed by a different author, making the analysis blind to sprouting location. Predicted sprouting location was restricted to the bottom 50% of pressure values (Fig. 5A). Within the lower pressure region, we identified all shear stress minima and all locations where vorticity equalled zero (Fig. 5B). Excluding locations with zero vorticity, the lowest minimum was chosen as the predicted sprout location (Fig. 5B). The actual sprout location was then assessed (Fig. 5C), blind to the predicted sprout location, based on the green images of vascular morphology and the location of x=0 from the pressure contour plots. In all embryos (n=4/4 embryos, P=0.001), we were able to correctly predict sprout location. Without VEGF, however, sprouts would not form at all. It is important to note that, in our system, fluid dynamics did not predict which avascular regions would sprout, but instead the flow predicted where along the vessel wall sprouting would occur. We did find that the sprouting avascular regions were larger on average than non-sprouting avascular regions (Fig. S2A, n=10 for each, P=0.012). Given the observed spread within the data (Fig. S2B), the size of the avascular region alone cannot predict which regions will sprout.

We frequently observed that sprouts would extend and retract before finally connecting to form a lumenised vessel (Fig. 6A, yellow arrows indicate direction of sprout movement), as has been observed in other model systems (Murakami et al., 2006). We plotted the elongation rate of the sprouts (Fig. 6B), measured by the displacement of the tip of the sprout between two frames, and compared this to changes in the three analysed hemodynamic parameters in our system. We observed that regression of the sprout occurred when the difference in pressure between the vessels became negative (Fig. 6C, n=3/3 embryos). We therefore plotted the elongation rate with respect to the pressure difference between the two vessels (Fig. 6D). Negative elongation rates represent sprout regression. We found a logarithmic relationship between the elongation rate and the pressure difference between the vessels (Fig. 6E). The non-linear regression gave an r² value of 0.82 (Fig. 6D), whereas the linear regression on the semi-ln plot gave an r² value of 0.87 (Fig. 6E, P<0.05). Thus, these results show that there is a significant association between the rate of elongation and the pressure difference between two vessels.

In some time-lapse experiments, we observed the direction of elongation changing dynamically in relation to pressure differences between two vessels. This occurred in two of the three arterial time-lapse recordings. In one example, taken from Embryo 2 of the arterial time-lapse recordings (Fig. 7), the sprout was extending at time point 210 min (Fig. 7A, yellow arrowhead). The middle panels show the contour maps for relative pressure in this region with the extending sprout labelled in black (Fig. 7B). At 210 min, the pressure difference between the vessels in the sprouting direction was positive (Fig. 7C). The flow pattern changed, resulting in a negative pressure difference in the direction of sprouting by 240 min. The sprout retracted slightly and then changed direction towards a point that resulted in a positive pressure difference (t=270 min, Fig. 7C). Sprouts did not necessarily extend along the maximum pressure difference. They would continue in the same direction as long as the pressure difference was positive but changed direction or regressed when the pressure difference was negative. As such, the pressure difference was permissive rather than instructive.

We show that flow dynamics have a significant effect on the sprouting location, the direction of sprouting and the rate of elongation in arterial vascular networks. We found that sprouts only form from a vessel at low pressure towards a vessel at higher pressure, such that there is a positive pressure difference at the sprouting location. Within the lower pressure region, sprouts formed at the lowest minimum in shear stress once locations where streams merge were excluded.

Our results do not, however, show that fluid dynamics control when a sprout will form, but rather where they are likely to form. Previous in vitro work on sprouting and shear stress had shown that shear stress could inhibit VEGF-induced sprouting (Song and Munn, 2011). Similarly, our results suggest that the formation of a sprout within an avascular region is initiated by non-hemodynamic parameters, but that the fluid dynamic parameters define the position of that sprout.

Our results support that shear stress is inhibitory to sprouting, as we and others have found (Chouinard-Pelletier et al., 2013; Song and Munn, 2011; Tressel et al., 2007). More recently, another group published that shear stress induced rather than inhibited sprouting (Galie et al., 2014). It is possible that both extremes are permissive to sprouting. Both low (<10 dyn/cm²) and pathologically high shear stress (>50 dyn/cm²) result in increased cell turnover by stimulating both apoptosis (Davies et al., 1986; Dolan et al., 2011) and proliferation (Davies et al., 1986; Metaxa et al., 2008; Sho et al., 2003). Physiological levels of shear stress, by contrast, induce endothelial cell quiescence. Our analysis indicates that the physiologically relevant signal during development is the presence of lower shear stress levels. Under pathologically high shear stress, we cannot exclude that shear stress could be pro-angiogenic.

Our previous results had indicated that shear stress only inhibited sprouting in veins but not in arteries (Chouinard-Pelletier et al., 2013). Similarly, others had found that steady shear stress inhibited sprouting in venous or capillary endothelial cell types (i.e. HUVECs and HMECs) but not in arterial endothelial cell types (BAECs) (Tressel et al., 2007). Our current results, however, indicate the opposite. Our previous results were based on an embryonic model in which we altered shear stress levels by altering the viscosity of the blood (Chouinard-Pelletier et al., 2013). Increased viscosity not only changes shear stress but would also create a more drastic pressure drop within the vascular network, such that pressure differences between vessels in the veins might be smaller. An alternative explanation is that increasing the shear stress levels throughout the embryo resulted in an arterialisation of the venous network. Though we cannot differentiate between these possibilities, it remains clear that the arterial and venous identity of the bed affects the relationship between shear stress and sprouting dynamics.

Our results also indicated that though locations where two flows merge create points of low shear stress and large pressure differences, sprouts do not initiate from these points. We investigated flow during the cardiac cycle to identify whether flow reversals or oscillations were present at these locations, but found that the flow remained unidirectional throughout the cardiac cycle. Whether vorticity affects endothelial cells has never been studied, but it seems improbable as vorticity is not a force but rather a measure of the rate of rotation of the fluid elements. The endothelium is a continuous layer connected by gap junctions and endothelial cells can sense gradients in shear stress (Dolan et al., 2011; Rouleau et al., 2010). Points at which there is a change in sign in the vorticity are surrounded by shear stress vectors on either side that are in opposite directions. As such, these points do not have a zero gradient in shear stress (i.e. where the magnitude and sign of shear stress is the same on both sides of the sprouting location) and do not represent true minima in shear stress.

We find that the pressure difference between vessels, and not the absolute value of the relative pressure, is proportional to elongation rate. It is improbable that the hydrostatic pressure itself is biologically active in our model. The pressure differences that we observed in these vessels are several orders of magnitude smaller than the levels of hydrostatic pressure that have been shown to be biologically active (Sato and Ohashi, 2005). The pressure difference between the vessels does, however, control the rate at which fluids exit a vessel and therefore regulates the interstitial flow patterns. Interstitial flow alone, without a VEGF gradient, can promote sprouting in an in vitro system (Song and Munn, 2011). Interstitial flow will modify the distribution of growth factors and metabolites in the interstitial tissue, and it was previously shown that a gradient of VEGF in the presence of interstitial flow was more effective at promoting sprouting than either alone (Song and Munn, 2011). As such, rather than sprout elongation being controlled by pressure directly, current in vitro evidence would indicate that the interstitial flow is the instructive parameter. We are currently developing computational models to study the distribution of angiogenic cues and the range of interstitial flow to further investigate this phenomenon.

All experimental methods were performed in accordance with the relevant guidelines on animal use. A detailed description of the methods was previously presented (Ghaffari, et al., 2015). Fertilised quail eggs (Cortunix japonica) were incubated at 37°C and ∼60% humidity until they reached a developmental stage of 12-14 somites or Hamburger Hamilton 11 (HH11). The embryos were injected intravascularly with labelling dyes using a picospritzer III micro-injector (General Valve Corporation). Alexa Fluor 488 acetylated low-density lipoprotein (AF488-AcLDL, Invitrogen, L-23380) was used to label endothelial cells and PEGylated polystyrene microspheres (0.5 µm diameter red, Invitrogen, Amino-PEG, Laysan Bio) were injected to image the flow dynamics. Embryo culture for time-lapse microscopy was done as previously described (Al-Roubaie et al., 2012). Embryos were imaged on an upright fluorescence microscope equipped with an Axiocam MRC for images of the vascular morphology (AF488-AcLDL) and a high-speed camera (Photron FASTCAM Ultima APX-RS) for images of the microsphere motion. Endothelial cells were imaged with a 5× objective lens every 15 min. Blood flow dynamics were imaged with a 10× objective lens at 250 fps for two full cardiac cycles (equivalent to 2 s of imaging time), once every 2 h.

The green image of the vasculature and the 500 red images of the microsphere motion were flattened into a single image to identify the walls of perfused vessels. A median filter was applied to the flattened image and then the image was thresholded to obtain a binary image used as the ‘image mask’ (for additional details, see Ghaffari et al., 2015). The image mask was applied to each individual frame of the microsphere motion. The centreline of each vessel was found by skeletonising the image mask. A Canny algorithm for edge detection was applied to the image mask to identify vessel walls. The local vessel diameter was obtained based on the vessel centreline and the vessel walls. Flow inlets and outlets were chosen far from branch points. Micro-PIV was used to calculate velocity waveforms based on the microsphere motion during the cardiac cycle at these inlets and outlets. Because the Reynolds number for the flow is very close to one and the Womersley number is less than one (Jones et al., 2004), entrance length effects are on the order of microns and a fully developed velocity profile was assumed for the µPIV analysis.

The vascular geometry obtained by image processing and the inlet/outlet velocities acquired by µPIV were imported into a computational fluid dynamic (CFD) solver. The governing equations inside the vessels (i.e. continuity and momentum equations for laminar, non-Newtonian, incompressible flow) were solved using a finite elements method (FEM) in COMSOL Multiphysics 4.4. The unsteady form of the momentum equations was used for flow inside the vessels to take into account the effects of blood flow pulsatility. The viscosity of blood was adjusted for embryonic stage, vessel diameter and shear rate as previously described (Ghaffari, et al., 2015). From N inlets and outlets, the flow rates at N−1 points were set as boundary conditions. For the last outlet, the fully developed boundary condition was applied at the extended outlet of the domain (dependent on Reynolds number), i.e. (∂/∂_x)=0. Mesh generation and solution was performed by COMSOL. The domain was meshed with extremely fine free triangular meshes. Two convergence criteria were used: first, a mass flux residual of less than 10⁻⁸ for each control volume; second, (|φ_i+1−φ_i|)/|φ_i+1|≤10⁻¹⁰ for all time steps where φ represents v_x or v_y (flow velocities in x and y directions) and i is the number of iterations.

For formation of sprouts on the lower pressure side, a chi-squared test was used to compare the probability of sprouts always forming on one side (n=6) with the probability of sprouts forming randomly (n=3 and n=3, for low and high pressure side, respectively). For sprout location at shear stress minima, the arc length was divided into six possible sprouting locations and a chi-squared test was used to evaluate the probably that sprouts would form in one location rather than randomly in all six locations. The value of six was chosen as a conservative value as it would equate to ∼30 µm of arc length. Sprouts form within microns of predicted shear stress minima. For the relationship between pressure difference and elongation rate, we tested for co-integration using a Dickey–Fuller test and established that the relationship between the variables was stationary (P<0.05).