Research articlesComputational simulations of the 4D micro-circulatory network in zebrafish tail amputation and regeneration
Abstract
Wall shear stress (WSS) contributes to the mechanotransduction underlying microvascular development and regeneration. Using computational fluid dynamics, we elucidated the interplay between WSS and vascular remodelling in a zebrafish model of tail amputation and regeneration. The transgenic Tg (fli1:eGFP; Gata1:ds-red) zebrafish line was used to track the three-dimensional fluorescently labelled vascular endothelium for post-image segmentation and reconstruction of the fluid domain. Particle image velocimetry was used to validate the blood flow. Following amputation to the dorsal aorta and posterior cardinal vein (PCV), vasoconstriction developed in the dorsal longitudinal anastomotic vessel (DLAV) along with increased WSS in the proximal segmental vessels (SVs) from amputation. Angiogenesis ensued at the tips of the amputated DLAV and PCV where WSS was minimal. At 2 days post amputation (dpa), vasodilation occurred in a pair of SVs proximal to amputation, followed by increased blood flow and WSS; however, in the SVs distal to amputation, WSS normalized to the baseline. At 3 dpa, the blood flow increased in the arterial SV proximal to amputation and through anastomosis with DLAV formed a loop with PCV. Thus, our in silico modelling revealed the interplay between WSS and microvascular adaptation to changes in WSS and blood flow to restore microcirculation following tail amputation.
1. Introduction
Biomechanical forces modulate vascular morphological adaptation during development [1,2]. While mechanotransduction mechanisms underlying arterial inflammatory responses and oxidative stress are well elucidated at the macro-scale in association with high Reynolds number flow [3], much less is known about the interplay between haemodynamic forces and vascular remodelling at the micro-scale in association with low Reynolds number flow [4,5]. Further, while endothelial dysfunction in the microvasculature is well recognized in myocardial ischaemia, diabetes and peripheral artery diseases (PADs) [6–8], investigation into vascular injury-mediated changes in haemodynamic forces and adaptive vascular remodelling in the microvasculature remained an experimental challenge.
Zebrafish (Danio rerio) embryos are transparent for light microscope imaging and genetically tractable for studying shear stress-mediated myocardial and vascular development [9]. We previously established four-dimensional (4D) light-sheet fluorescent microscopy (LSFM) to study the development of cardiac trabeculation during zebrafish cardiac morphogenesis [10,11]. Combining LSFM and computational fluid dynamics (CFD), we have provided biomechanical insights into wall shear stress (WSS)-mediated endocardial trabeculation and outflow tract valvulogenesis in embryonic zebrafish [12–14]. Thus, we sought to use the transgenic zebrafish Tg(fli1:eGFP; Gata1:ds-red) line, in which enhanced green fluorescent protein-labelled endothelium allows for segmentation of the three-dimensional (3D) vasculature, and ds-red-labelled blood cells (mainly erythrocytes) for CFD reconstruction and particle image velocimetry (PIV).
WSS intimately regulates the mechanotransduction underlying vascular development, repair and homeostasis [15,16]. The spatial and temporal variations in fluid shear stress in the segmental vessels (SVs) of zebrafish embryos modulate the transformation of arterial to venous vessels, leading to a balanced distribution of arterial and venous SVs during vascular development [17]. In this context, the zebrafish microvascular network can be considered as a circuit analogy model consisting of the dorsal aorta (DA), SVs and the posterior cardinal vein (PCV). This model allows the microvascular response to blood cell flux partitioning in each SV to be simulated, leading to uniform flow distribution [18,19]. However, this zero-dimensional resistance modelling neglects bidirectional flow in the dorsal longitudinal anastomotic vessel (DLAV) network, vascular remodelling, including constriction and dilation, and angiogenesis of the DLAV and PCV to form a new loop following tail amputation. Zero-dimensional modelling approaches also overlook the local variations in vessel diameter in association with the haemodynamic resistance in the microvasculature. For this reason, our CFD modelling sought to account for haemodynamic resistance by using the detailed morphology of the tail and performing a full simulation for accurate estimation of WSS and luminal interaction.
We hereby demonstrate a 3D CFD model encompassing four time-lapse stages of vessel regeneration following tail amputation, and subsequent loop formation between the DA and PCV. We employed eGFP-labelled vascular endothelium by sub-micron confocal imaging and ds-Red-labelled blood cells to reconstruct the computational model and to simulate the interplay between haemodynamic changes and vascular morphological remodelling. We discovered a coordinated vascular adaptation to tail amputation that disconnected the DA from the PCV, followed by DLAV vasoconstriction at 1 day post amputation (dpa), and SV dilation at 2 dpa, coupled with the rise and fall of WSS and the haemodynamic changes in flow rates. At 3 dpa, the flow rate in the arterial SV proximal to amputation continued to rise and the arterial SV merged with DLAV to form a new loop with PCV. The microvascular remodelling in the zebrafish tail injury model comprises arteriogenesis in SVs proximal to the amputation site followed by sprouting angiogenesis to replace the main vascular network (DA → PCV replaced by DA → arterial SV → regenerated vessel → PCV). Thus, elucidating the interplay between shear forces and remodelling provides new biomechanical insights into the micro-scale haemodynamics, with translational implications in similar vascular injury and repair mechanisms in patients with diabetes or PAD undergoing revascularization [6,20–22].
2. Material and methods
Zebrafish (Danio rerio) has the capacity to regenerate following distal tail amputation, providing a genetically tractable model to study changes in microcirculation in response to vascular development, injury and regeneration [23–25]. We performed tail amputation on embryonic zebrafish to assess the dynamic changes in WSS and flow rates in response to microvascular adaptation and remodelling proximal to the amputation site. The investigation of haemodynamic changes and morphological remodelling was established by imaging the transgenic Tg(fli1:eGFP; Gata1:ds-red) zebrafish embryos from 4 to 7 days post fertilization (dpf) (figure 1a–c), followed by post-imaging segmentation (figure 1d,e) and 3D blood flow simulation (figure 1g,h). The enhanced green fluorescent protein (eGFP)-labelled endothelial cells expressed vascular endothelial growth factor (VEGF) receptor, denoted as fli1, and the ds-red-labelled blood cells (mostly erythrocytes and other blood lineages) allowed for tracking their movement. Thus, sequential analyses from stitching the 3D confocal imaging to computational modelling and PIV were performed to validate the haemodynamic changes proximal to the amputation site. Microvascular constriction and dilation and new vessel formation (angiogenesis) in response to tail amputation were captured and implemented on the model to estimate haemodynamic variations (figure 1).
Figure 1. A pipeline of post-imaging segmentation, model creation, CFD simulation, validation and analysis. (a) Confocal microscopy was used to image the vascular system from a transgenic Tg (fli:eGFP) zebrafish embryo at 5 dpf. The vasculature was visualized via the enhanced green fluorescent protein (eGFP)-labelled endothelial cells that express VEGF receptor, fli1. (b) Stitching of the 3D confocal images provided a reconstruction of the vasculature for post-imaging segmentation. (c) The ds-red-labelled blood cells (gata1:ds-red) in the zebrafish tail microvasculature were monitored at two different time points. The motion of ds-red-labelled blood cells allowed for quantitative assessment of the blood flow at t = to and t = to + dt in the microvascular network. Following segmentation, the 3D computational domain was reconstructed (d) prior to and (e) after tail amputation. (f) PIV was performed to track the individually ds-red-labelled blood cells in the SV from t = to to t = to + dt. (g) SimVascular open-source software was employed with boundary conditions derived from PIV. (h,i) CFD analyses revealed WSS prior to (h) and after (i) tail amputation. The magnitude of velocity in SVs was validated. Scale bar: 100 µm.
2.1. Transgenic zebrafish model
The transgenic Tg (fli1:eGFP; Gata1:ds-red) embryos were rinsed with fresh standard E3 medium and cultivated at 28.5°C for 4 days after fertilization. Standard E3 medium was supplemented with 0.05% methylene blue (Sigma Aldrich, MO) to prevent fungal outbreak and 0.003% phenylthiourea (PTU; Sigma Aldrich, MO) to reduce melanogenesis. At 4 dpf, embryos were randomly selected for immobilization with neutralized tricaine (Sigma Aldrich, MO) to undergo tail amputation. Zebrafish responses to anaesthesia were monitored by the movement of the caudal fin. The tail was amputated at approximately 100 µm from the distal end. Following amputation, embryos were returned to fresh E3 medium and maintained for 3 days.
2.2. Post-amputation imaging of blood flow and microvasculature
The images of blood cells (Gata1:ds-Red) transport were captured using an inverted microscope (Olympus, IX70) and digital CCD camera (QIclick, Teledyne Qimaging, Canada). We performed time-lapsed imaging from 4 dpf to 7 dpf or from 0 dpa to 3 dpa for both the amputated and control fish (amputation = 7; control = 5). The acquired images revealed the motion of erythrocytes in the microvasculature throughout multiple cardiac cycles.
The transgenic zebrafish Tg (fli1:eGFP; Gata1:ds-red) line was employed to capture GFP-labelled endothelium lining of the microvasculature undergoing remodelling following amputation from 4 to 7 dpf. We used a confocal imaging system (Leica TCS-SP8-SMD) to provide isotropic sub-micron resolution in three dimensions (approx. 500 nm in-plane and in the z-direction) for segmentation of the 3D microvasculature. To reach this resolution, we used the high-magnification (64× objective lens) images at 4–7 dpf, and we used 5 dpf for segmentation of the tail microvascular network.
We used our in-house light-sheet fluorescent microscope for the time-lapse imaging of zebrafish tails to assess the vascular remodelling in response to amputation [26]. Our system used a continuous-wave laser (473 nm; Laserglow Technologies, Canada) and a dry objective (Plan Fluor 4×/0.13; Nikon) to generate a planar illumination. A detection module including a water dipping objective (Fluor 20×/0.5w; Nikon) and a scientific CMOS camera (ORCA-Flash4.0; Hamamatsu, Japan) was positioned perpendicular to the illumination. Our system provides a resolution of 2 µm, which was adequate for the measurement of microvascular diameter variation after amputation [26].
2.3. Segmentation of the microvascular network
We used the 3D GFP-labelled images of endothelial lining to reconstruct the computational domain in SimVascular (figure 1b) [27]. The pathline of each vessel was produced by connecting the points in the centre of the lumen along the vessel in the 3D image stack. For each vessel, the cross-section was calculated by interpolation of the manually specified points on the wall/lumen interface (one dimension → two dimensions). Linking the cross-section areas along the vessel path with a spline function, we reconstructed each vessel in the network (two dimensions → three dimension). Overall, we generated the vascular network for 21 SVs and 20 DLAVs, DAs, PCVs and capillary venous plexuses (figure 1d,e).
We built the post-amputation model by morphing the geometrical changes to the pre-amputation model (figure 2a–c). At 1 dpa, we removed eight SVs (six SVs for amputation and two SVs for vasoconstriction) from the zebrafish tail microvascular network distal to the heart to model the microvascular response to amputation followed by vasoconstriction at the DLAV (1 dpa). At 2 dpa, two SVs were added back to the microcirculation network where vasodilation was implemented (figure 2d,e). Confocal and light-sheet images of eight embryonic zebrafish post amputation were used for measurement of vasodilation magnitude. The measurement was performed at more than 10 locations along each SV perpendicular to the SV profile in the plane of the SV section (electronic supplementary material, figure S4). The diameter values are averaged along the vessel and then averaged among eight zebrafish. We employed the time-lapse LSFM images to measure the average length and diameter of the regenerated vessel. The arterial SVs proximal to amputation were connected to the PCV to model loop formation at 3 dpa (figure 2c).
Figure 2. Zebrafish embryonic cardiovascular circulatory system in response to tail amputation and regeneration. (a) Arterial circulation is denoted in red and venous in blue. Two regions were magnified to reveal the vascular network; namely, top view and tail amputation. (b) The top view reveals that the DA originates from the ventricle, and the PCV returns to the atrium. Arterial SVs connect with the DLAV, which drains to the PCV via venous SVs. (c) The tail amputation site (side view) illustrates SVs undergoing vasoconstriction, vasodilation and loop formation proximal to amputation at 1, 2 and 3 dpa. The proximal SVs were numbered from 1 to 6. (d) The changes in proximal SV diameters were compared from 0 dpa to 3 dpa. Both SV1 and SV2 underwent significant changes in diameter due to vasoconstriction at 1 dpa and vasodilation at 2 dpa (p < 0.05 versus control, n = 5). The changes at 3 dpa were statistically insignificant. (e) Changes in diameters were significant to the proximal SV1 and SV2 as compared with the distal SV3–6. n.s., not significant; DA, dorsal aorta; SV, segmental vessel; DLAV, dorsal longitudinal anastomotic vessel; and PCV, posterior cardinal vein; dpa, days post amputation.
2.4. Particle image velocimetry for the boundary conditions and validation
We developed and implemented our PIV code on ds-Red-labelled blood cells. To determine the boundary conditions for the inlet (DA) and to validate the results, we measured the velocity in the DA and individual SVs. The velocity was measured over 20 cardiac cycles in each individual vessel and the time-averaged value was used for the inlet boundary condition at DA and for validation at SVs that were visible in blood flow images. Our custom-written PIV algorithm was able to correlate with the intensity distribution of ds-Red-labelled erythrocytes in different time frames. A two-dimensional continuous polynomial function was fitted to the correlation distribution, and the position of the maximum correlation was determined from the polynomial function between two time frames. At 20 frames per second, the computed erythrocyte position allowed for estimating the velocity in the DA and SVs in which the blood cells were visible (figure 1c,f). The velocity profile in SV pairs at the same distal location was filtered based on the curvature of SVs and blood flow direction (arterial versus venous) in them to distinguish the arterial and venous SVs at different depths.
2.5. Three-dimensional computational modelling
A stabilized 3D finite-element method was applied for CFD analysis to simulate the blood flow in the microvascular network using SimVascular's svFSI solver [27,28]. The Navier–Stokes equations were solved in svFSI with an incompressible flow assumption to quantify haemodynamic changes in response to tail amputation and vascular regeneration from 0 to 3 dpa (4–7 dpf) (figure 1g). In the setting of low Reynolds number (, where D is defined as the vessel diameter) in the microvascular network, viscous forces dominate over inertial forces, and the diameters and lengths of the vessels governed the pressure distribution and flow rates. Owing to the low Reynolds number ( , where D is defined as the erythrocyte diameter), the pathway of the ds-Red-labelled erythrocytes followed the streamline of blood flow. Unlike anucleated red blood cells in the mammalian circulation, zebrafish erythrocytes contain nuclei and are, thus, less deformable.
Despite the pulsatile nature of blood flow in zebrafish microcirculation, the Womersley number (ratio of transient inertial force to viscous force: ) was small, and the time-dependent terms in the Navier–Stokes equations were negligible in our continuum modelling [27]. The flow rate (, where Q is the flow rate, V is the velocity and A is the cross-sectional area) in each SV was calculated by integrating the velocity in a cross-section. The WSS ( where μ is viscosity and r is the radius) was calculated and averaged along the vessel wall in each SV.
We calculated the inlet velocity into the DA during a cardiac cycle using the PIV code, and the time-averaged magnitude of Vave = 238 µm s−1 was applied as the inlet boundary condition. Given the small Womersley number (α ≪ 1), a fully developed parabolic flow was considered in the inlet of the DA (red arrow in figure 1f). A zero-pressure boundary condition was applied at the outlet of the PCV (blue arrow in figure 1f). A non-slip boundary condition was assumed on the vessel walls.
To incorporate non-Newtonian effects of flow in the microcirculation, we used a Carreau–Yasuda viscosity model. The parameters of this model were extracted from our previous study on non-Newtonian behaviour of zebrafish blood [29] in a power law model,
3. Results
3.1. Segmentation of the three-dimensional microvasculature
Following 3D stitching of the confocal images from the intact and amputated zebrafish tail, we performed segmentation of the microvascular network (figure 3a,b). Next, the 3D computational domain was reconstructed to simulate haemodynamic changes in response to tail amputation and vascular regeneration. The box indicates the region of interest for vascular segmentation and simulation (figure 3a). A representative 3D computational domain was acquired from the segmentation process (figure 3c,d). The pressure drop across the venous plexus was deemed to be negligible owing to both the short vessel length and the low haemodynamic resistance. These considerations allowed for connecting the venous SV with the PCV for segmentation through the largest venous plexus.
Figure 3. (a) A transgenic Tg(fli1:GFP) zebrafish line was used to visualize the 3D vascular network. The dashed rectangle (yellow) indicates the tail region from which vascular network was reconstructed for CFD simulation. (b) Magnification of the tail region depicts the flow direction in the DA in relation to the PCV. (c) Representative segmentation of the 3D vasculature recapitulates the DA, SV, DLAV and PCV. (d) Magnification of the yellow rectangle reveals a cross-section of the tail region where the vascular connectivity occurs between the DLAV and PCV and between the DA and SVs. Note that the arterial SVs branch off from the DA, whereas venous SVs drain into the PCV. CFD: computational fluid dynamics; DA, dorsal aorta; SV, segmental vessel; DLAV, dorsal longitudinal anastomotic vessel; and PCV, posterior cardinal vein. Scale bars, 50 μm.
While DLAVs are commonly assumed to be a single, straight vessel for unidirectional flow, bidirectional flow occurs in the DLAVs in parallel. This is demonstrated in ds-Red-labelled blood cell images (electronic supplementary material, video S1) and the track of particles on CFD-derived streamlines (electronic supplementary material, video S2). For this reason, a parallel pattern was proposed for the DLAVs according to their interconnection in the reconstructed vascular network (figure 2b; electronic supplementary material, figure S1). Thus, we demonstrate a pair of DLAVs in parallel, intersecting with the arterial and venous SVs (in blue and red) to regulate the blood flow between SVs.
3.2. Vascular remodelling proximal to the tail amputation site
Confocal and light-sheet images revealed remodelling of the vessels and changes in blood flow from 0 to 3 dpa (figure 2). At 0–1 dpa, DLAVs underwent vasoconstriction, accompanied by a reduction in blood flow in the SVs proximal to the amputation site while both the DLAV and PCV started to form a new connection towards the tail region (electronic supplementary material, figure S2 and videos S3 and S4) [31]. At 0–1 dpa, SVs proximal to the amputation site (SV a + 1, a + 2, where a denotes the number of SVs that were removed after amputation) underwent statistically significant vasodilation (SV a + 1 and SV a + 2, p < 0.0001, n = 61), whereas the diameter variations in the SVs distal from the amputation site (SV number > a + 3) remained statistically unchanged (SV a + 3 to a + 6, p > 0.05, n = 61) (figure 2d,e; electronic supplementary material, video S4). At 2 dpa, the DLAV proximal to SV a + 1 started to form a new loop with PCV (figure 2c; electronic supplementary material, videos S5 and S6). At 3 dpa, vasodilation in SV a + 1 and a + 2 attenuated in conjunction with the complete loop formation between the DLAV and PCV, allowing for vascular regeneration to restore blood flow in the tail (figure 2c; electronic supplementary material, videos S7 and S8). Thus, we unveil micro-circular remodelling proximal to the tail amputation site, where the diameters of the DLAV reduced at 1 dpa, followed by increased diameters in arterial SV a + 1 and a + 2, with concomitant angiogenesis in the DLAV and PCV to form a new loop at 2 dpa, followed by normalization of SV diameter when a new loop connected the arterial and venous circulation.
3.3. Particle image velocimetry to validate computational fluid dynamics
To validate our CFD model for low Reynolds number flow in the tail region, we performed PIV to obtain the mean velocity of blood flow in the SVs from proximal to distal to the amputation site. CFD simulation revealed velocity ranging from 130 to 550 µm s−1, which was comparable to that obtained by PIV (figure 4e). The CFD results deviated from those of PIV by approximately 8%, and the velocity magnitude of the most distal SV (number 15) was significantly higher than that of the proximal SVs, indicating a good agreement with the experimental data.
Figure 4. CFD simulates changes in haemodynamic responses to tail injury and regeneration. Velocity streamlines in the computational domain were colour-coded by the velocity magnitude (a) at 4 dpf prior to tail amputation or 0 dpa, (b) 1 dpa, (c) 2 dpa and (d) 3 dpa. Arrow indicates vascular repair with a loop formation. (e) Variations in velocity from the proximal to distal SV were compared between CFD and PIV prior to tail amputation. (f) The dynamic changes in flow rates developed from the proximal to distal SVs before (4 dpf or 0 dpa) and after tail amputation (5–7 dpf or 1–3 dpa). Note that the flow rate in SV 7 increased at 3 dpa when DLAV formed a loop with PCV.
3.4. Changes in the haemodynamic parameters prior to tail amputation
Results from the 3D CFD simulations included WSS, pressure gradients, velocity profiles and flow rates in the individual SVs. Our model verified the unique pattern of blood flow that recirculates between two SVs in parallel (electronic supplementary material, videos S1 and S2). The blood velocity and flow rate in the DA decreased towards the tail region where arterial blood flow drains into the SVs and DLAVs. Thus, this compartmental increase in haemodynamic resistance in both the SVs and DLAVs towards the tail region results in a decrease in flow rate in the DA.
Our 3D simulations further demonstrated the bidirectional flow in the DLAVs in parallel (see figure 2b; electronic supplementary material, video S2). The flow direction in each DLAV was regulated by the pressure difference between the arterial and venous SVs corresponding to the diameter distribution in SVs under Stokes flow . As a result, the combination of DLAVs in parallel and SVs of various diameters modulated the flow rates and WSS (figure 5d). The flow rate in each SV was primarily regulated by the changes in diameters of the arterial and venous SVs and the interconnecting DLAV (figure 4d,f). Taken together, our in silico model predicted that small changes in the arterial SV diameter and its branching vessels regulate the direction of blood flow in the microvascular network.
Figure 5. CFD-simulated changes in WSS in response to tail amputation and regeneration. (a) Relatively high WSS developed in the distal arterial SV at 4 dpf prior to tail amputation. (b) An increase in WSS occurred in the proximal SV following tail amputation at 1–2 dpa. (c) An increase in WSS developed in the DLAV and arterial SV proximal to the amputation site, along with a reduction in WSS in the rest of the SVs at 2 dpa. (d) A reduction in WSS developed in distal SVs in response to a loop formation connecting DLAV with PCV at 3 dpa. (e) The mean WSSs from the distal SVs were compared prior to and post tail amputation, followed by loop formation at 3 dpa.
3.5. Changes in wall shear stress in response to tail amputation and regeneration
We compared changes in flow rates and WSS in the vascular network from 4 dpf (pre-amputation) to 7 dpf (3 dpa—loop formation) (figures 4 and 5). Next, we compared the changes in the area-averaged WSS (mean WSS) and flow rates in SVs following tail amputation (figure 5b,c). At 0 to 1 dpa, we used the transgenic Tg (fli1:eGFP) zebrafish line where VEGF receptors expressed in the vascular endothelial cells were encoded with the eGFP. We were able to visualize the DLAVs proximal to the amputation site undergoing vasoconstriction, restricting blood flow, along with a reduction in WSS in SVs 7, 8 (electronic supplementary material, figure S2 and video S2). Our simulation further showed that WSS increased following tail amputation in the rest of the SVs between 0 and 1 dpa (figure 5a–e). This increase was preferentially higher in the arterial SVs (i.e. SVs 9, 12) than in the venous SVs (i.e. SVs 10, 11) (figure 5e). The differentially higher levels of increase in flow rate at arterial SVs 9, 12 (compared with venous SVs 10, 11) is due to the diameter difference (DSV9,12< DSV10,11). In other words, the same magnitude of increase in flow rate would lead to a higher increase of WSS in those SVs with smaller diameter. The sum of flow rate variation in the arterial SVs versus venous SVs is equal owing to the mass conservation law at DLAVs. At 0–1 dpa, the eGFP-labelled images revealed that DLAV and PCV formed 75% of the new vessel length while proximal SVs (SVs 7, 8) developed a relatively low WSS owing to DLAV vasoconstriction and the average WSS in the rest of the SVs was elevated compared with unamputated fish (figure 5a–e).
Simulation of blood flow at 2 dpa predicted that vasodilation of the proximal SVs (SVs 7, 8) and subsequent new loop connecting DLAV with PCV resulted in a proximal increase in flow rates and normalization in WSS in the vascular network to the range prior to amputation (figure 5c,e; electronic supplementary material, video S5). However, WSS and flow rates in the distal SVs decreased at 2 dpa to forestall further SV dilation (figure 5c,e). Our simulation also predicted that vasodilation increased the flow rate proximal to the amputation site (SV, 7, 8), while the reduction in WSS further attenuated SV dilation (see Material and methods). Thus, our model demonstrates that SV dilation mediated a reduction in WSS (figures 2f and 5e).
The CFD model captured the loop formation of the regenerated vessel occurring distal to the arterial SV proximal to the amputation site (SV 7). This SV supplied the arterial flow to DLAV, which underwent angiogenesis to form a new loop with PCV (figures 4d and 5d). This flow path (DA → SV 7 → DLAV → regenerated vessel → PCV) decreased the blood supply to the rest of the SVs and DLAVs, leading to the DA → PCV loop formation, and an increase in flow rates and WSS in the arterial SV proximal to amputation (figures 4d,f, and 5d,e). However, the blood flow to the venous SV proximal to the amputation site was reduced by the new loop between DLAV and PCV, resulting in a reduced WSS (figures 4d,f and 5d,e). Overall, CFD simulation demonstrated the changes in flow rate and WSS in response to the variation in SV diameters and flow rates and indicated microcirculation following tail amputation and regeneration.
3.6. The non-Newtonian effects on the flow rates and wall shear stress
From the power law model, we assessed the effects of the non-Newtonian index on the microcirculation. We normalized the results according to the reported value for the non-Newtonian index (n = 1.09). The normalized mean WSSs in SVs were compared with different non-Newtonian indices (figure 6a,c). Our parametric study indicated a nonlinear effect of the non-Newtonian index on WSS, with n = 1.09 indicating a maximum WSS condition (figure 6c). Similarly, the flow rates in SVs exhibited a comparable trend (figure 6b), and the overall circulation of blood was independent of the non-Newtonian behaviour.
Figure 6. The effects of the non-Newtonian index (n) on WSS contours from 1 to n = 1.09. (a) The non-dimensionalization was performed using the experimental index at n = 1.09. Variations in WSS and flow rates were nonlinear. The highest value for WSS was obtained using the standard value n = 1.09. (b) Dimensionless WSS and (c) flow rates were compared from proximal to distal SVs.
We compared WSS contours for values of n from 1 to 1.09 (figure 6a). We indicated that flow rate is insensitive to the non-Newtonian model indices.
4. Discussion
Fluid shear stress imparts both metabolic and mechanical effects on vascular endothelial function [3]. During development, haemodynamic forces such as shear stress are intimately linked with cardiac morphogenesis [32]. We have demonstrated that peristaltic contraction of the embryonic heart tube produces time-varying shear stress (∂τ/∂t) and pressure gradients (∇P) across the atrioventricular canal in a transgenic zebrafish Tg (fli-1:eGFP; Gata1:ds-red) model of cardiovascular development [1]. The advent of the zebrafish genetic system has enabled the application of the fli1 promoter to drive expression of eGFP in all vasculature throughout embryogenesis [33]; thereby, allowing for 3D visualization for CFD simulation. By integrating 3D reconstruction of confocal images with haemodynamic simulation, we unravelled the changes in WSS and microvascular adaptation in response to tail amputation and vascular regeneration. Tail amputation resulted in disconnection between the DA and PCV, which triggered a well-coordinated DLAV vasoconstriction, accompanied by the initial decrease in WSS and flow rate in the SVs proximal to the amputation site. Simultaneously, an increase in WSS and flow rate in the SVs distal to the amputation site developed as a result of bypass of the blood from the DA to PCV during 0–1 dpa. WSS normalized to baseline in the SVs during vasodilation at 2 dpa when the DLAV and PCV formed a new loop at 3 dpa, whereas flow rates in arterial SV merged with those in the DLAV and continued to rise. The size and velocity scales and therefore Reynolds number in our model are comparable with the microcirculation in other microvascular beds, including rodents and human (we used tricaine to anaesthetize the fish for imaging, which theoretically reduces the blood flow) [34–36]. Thus, our data provide new haemodynamic insights into microvascular network injury and repair at low Reynolds numbers.
Haemodynamic shear forces are well recognized to modulate endothelial homeostasis [37,38], migration [39], vascular development and regeneration [40]. At a steady state, shear stress (τ) is characterized as the dynamic viscosity (μ) of a fluid multiplied by its shear rate , defined as the tangential velocity gradient [3]: . The zebrafish system allows for genetic manipulations of viscosity to alter shear stress [12,25,41]. We previously reported that injection of Gata-1a oligonucleotide morpholino (MO) to inhibit erythropoiesis resulted in a decrease in viscosity-dependent shear stress and subsequent impaired vascular regeneration following tail amputation [42]. After TNNT-2 MO injection to inhibit myocardial contraction and blood flow, zebrafish embryos also failed to develop tail regeneration at 3 dpa, whereas the injection of erythropoietin (Epo) mRNA activates erythrocytosis to increase viscosity, resulting in restored vascular regeneration [42]. Thus, the rise and fall of WSS in adaptation to vascular remodelling in the SV and DLAV proximal to the amputation site occurred during the formation of a new loop connecting between the DLAV and PCV.
Our computational modelling was developed from the dynamic adaptation of vessel morphology to tail injury and vascular regeneration in the micro-scale under low Reynolds numbers. Choi et al. [19] studied the variation in WSS in an embryonic zebrafish tail using a circuit analogy model, and they demonstrated changes in WSS in the microvascular network by partially occluding the blood flow in SVs. Chang et al. [18] further proposed to incorporate the haemodynamic resistance of individual blood cells in the SVs towards the tail region. However, zero-dimensional resistance modelling does not capture vascular remodelling-associated changes in flow rates and WSS. Given the complex geometrical features of the zebrafish microcirculation, the variations in vessel diameter are pivotal for haemodynamics. For this reason, our CFD model addressed vasoconstriction and vasodilation, and vascular regeneration in the microvascular network, demonstrating a spatio-temporal-dependent modelling of vascular injury and repair in vivo.
Spatial and temporal variations in WSS modulate vascular endothelial function and homeostasis. WSS in the microcirculation regulates the canonical Wnt/β-catenin signalling pathway for vascular development [25,43], and WSS-activated Notch signalling modulates endothelial loop formation and endocardial trabeculation [12,23,42]. We have previously developed light-sheet imaging methods to study the development of trabeculation for zebrafish ventricles [10,11]. Combining light-sheet microscopy and CFD, we assessed the role of WSS in the formation of trabeculation and outflow tract valvulogenesis in embryonic zebrafish [12–14]. In this CFD modelling of vascular injury and repair, we further demonstrate the integration of simulation with the genetic zebrafish system to uncover microvascular mechanics.
Blood effective viscosity provides the rheological basis to compute WSS underlying cardiovascular morphogenesis. During development, blood viscosity regulates fluid shear stress to impart mechano-signal transduction to initiate cardiac trabeculation and outflow tract formation [12,32]. However, the minute amount of zebrafish blood renders measurements of blood viscosity experimentally challenging [44]. We previously reported the capillary pressure-driven principle to develop microfluidic channels for measuring the small-scale blood viscosity [29], allowing a power law non-Newtonian model to be established to predict the transient rise and fall of viscosity during embryonic development. In this context, we further analysed the non-Newtonian effects on the flow rates and WSS in SVs by performing a parametric study based on the power law index (n) (electronic supplementary material, figure S1). From the power law model, we normalized the results according to the reported value for the non-Newtonian index (n = 1.09) [29]. Our results indicate that, despite local variations in WSS, the blood flow in embryonic zebrafish tail was not significantly altered by variation in the non-Newtonian index. Overall, we unravelled and quantified the vascular remodelling and WSS variation in the face of tail amputation. Our CFD simulation was corroborated with PIV. The WSS-mediated arteriogenesis in SVs proximal to the amputation site is followed by sprouting angiogenesis to replace the main vascular pathway; thus providing the clinical translation to microvascular remodelling and revascularization in PADs that are associated with diabetes [6,20–22].
Ethics
The animal study was reviewed and approved by the University of California, Los Angeles (UCLA), under animal welfare assurance no. A3196-01.
Data accessibility
The data will be provided upon contact with the corresponding author.
Authors' contributions
M.R.: conceptualization, data curation, formal analysis, investigation, methodology, validation, visualization, writing—original draft and writing—review and editing; K.I.B.: conceptualization, data curation, formal analysis, methodology, visualization and writing—review and editing; Z.W.: data curation, formal analysis and writing—review and editing; S.C.: conceptualization, methodology, writing—review and editing; S.S.: writing—review and editing; A.L.: writing—review and editing; R.O.: writing—review and editing; V.V.: conceptualization, methodology, software, writing—review and editing; Y.D.: data curation, methodology and writing—review and editing; A.L.M.: methodology, software, supervision and writing—review and editing; T.K.H.: conceptualization, funding acquisition, methodology, project administration, resources, supervision, writing—original draft and writing—review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed herein.
Competing interests
We declare we have no competing interests.
Funding
The present work was funded by the National Institutes of Health to T.K.H. (grant nos. R01HL083015, R01HL111437, R01HL129727, R01HL118650 and VA MERIT AWARD I01 BX004356).
