A cell-and-plasma numerical model reveals hemodynamic stress and flow adaptation in zebrafish microvessels after morphological alteration

The development of a functional cardiovascular system ensures a sustainable oxygen, nutrient and hormone delivery system for successful embryonic development and homeostasis in adulthood. While early vessels are formed by biochemical signaling and genetic programming, the onset of blood flow provides mechanical cues that participate in vascular remodeling of the embryonic vascular system. The zebrafish is a prolific animal model for studying the quantitative relationship between blood flow and vascular morphogenesis due to a combination of favorable factors including blood flow visualization in optically transparent larvae. In this study, we have developed a cell-and-plasma blood transport model using computational fluid dynamics (CFD) to understand how red blood cell (RBC) partitioning affect lumen wall shear stress (WSS) and blood pressure in zebrafish trunk blood vascular networks with altered rheology and morphology. By performing live imaging of embryos with reduced hematocrit, we discovered that cardiac output and caudal artery flow rates were maintained. These adaptation trends were recapitulated in our CFD models, which showed reduction in network WSS via viscosity reduction in the caudal artery/vein and via pressure gradient weakening in the intersegmental vessels (ISVs). Embryos with experimentally reduced lumen diameter showed reduced cardiac output and caudal artery flow rate. Factoring in this trend into our CFD models, simulations highlighted that lumen diameter reduction increased vessel WSS but this increase was mitigated by flow reduction due to the adaptive network pressure gradient weakening. Additionally, hypothetical network CFD models with different vessel lumen diameter distribution characteristics indicated the significance of axial variation in lumen diameter and cross-sectional shape for establishing physiological WSS gradients along ISVs. In summary, our work demonstrates how both experiment-driven and hypothetical CFD modeling can be employed for the study of blood flow physiology during vascular remodeling.

Section B: Experiment flow data of trunk vessels from gata1 and control morpholino injected embryos at 2 dpf Fig A .RBC velocity data in the caudal artery CA (i), caudal vein CV (ii) and intersegmental vessels ISVs (iii) in trunk region pertaining to the simulation model domain (posterior to the yolk extension) for gata1 morphants at 2dpf and their corresponding RBC counts in the ISVs (iv).RBC velocity data in the CA (v), CV (vi) and ISVs (vii) in trunk region posterior to the yolk extension for control morpholino injected embryos at 2dpf and the corresponding RBC counts into ISVs (viii).Orange grid lines in velocity plots (i to iii, v to vii) indicate the time positions of systolic peaks in the CA velocity and show the phase lag between CA peaks and CV peaks.Red-dashed box in (v to viii) highlight Control MO #3 zebrafish, which was used as the normalization reference for calculating the relative CA hematocrit in both gata1 MO and control MO groups.

Section F: The coarse-grained spectrin model (CGSM) of RBC mechanics
The zebrafish RBC is modelled as a flattened ellipsoid capsule with an RBC membrane shell enveloping a cytosolic interior and a nucleus.The membrane is represented by a hexagonal mesh where edges represent the cytoskeleton (CSK) of spectrin chains joined to actin filaments [2,3] and the triangle surface elements represent the plasma membrane (PM).Likewise, the nucleus is defined as a three-dimensional meshwork where enjoining edges act as compressible springs for the prescription of nuclear deformation.Tension in the RBC membrane arising from deformation consists of the elastic component and the viscous component:  , =  , +  , (1) where  , is the viscous force and  , is the elastic force at membrane vertex .The elastic force may be obtained by differentiating the local elastic energy ( , ) with respect to the nodal displacement (  ) produced by the deformation: The elastic constitution of the RBC membrane can be decomposed to the individual mechanical aspects of the PM and CSK:

Shear elasticity of the CSK and stiffness of the RBC nucleus
Shearing properties of the CSK were modelled using the worm-like chain (WLC) extension model and the power-law (POW) compression model (Fedosov et al., 2010a): where   is the mesh edge length between vertices  and ,   is the maximum allowable edge length in the WLC model and  is the persistence length.   is the energy per unit mass which may be considered as a membrane constant for an isothermal consideration.  is the spring constant for the POW compression spring model.The resting shear modulus ( 0 ) of the WLC-POW network can be given by the following relation [4]: where  ,0 is the mesh edge length for the RBC at rest.
Using the WLC-POW springs we are able to exhibit the strain hardening behavior of RBC membranes (see Fig. 10A).We employed the same WLC-POW spring elements for the discretized network of springs that model the stiffness of the RBC nucleus.

Bending elasticity
The bending energy in the mesh network was calculated using the angle   subtended by the normals ( ) of the two adjoining triangular mesh elements along their shared mesh edge ij: where   is the bending modulus,  ,0 is the local spontaneous angle that a free sheet of spectrin polymers will adopt in their unstressed state.Since we assume the RBC rest state to be the reference energy state in our model,  ,0 is the angle between triangular mesh elements for the RBC at rest.

Area elasticity
The PM has been reported to be highly incompressible and area dilations exceeding 2-4 % produce immediate lysis [5,6].Incompressibility of the PM was prescribed by the area penalty model: where the area deformation energy is a combination of both global and local surface area changes.  is the current area of the local triangular mesh element and  ,0 is the area of the element when the RBC is at rest.  is the current surface area of the RBC membrane and  ,0 is the membrane surface area of the RBC at rest.  represents the global area compressibility coefficient while   is the local compressibility coefficient.The compressibility coefficients may be related to the area compressibility modulus by the following relation [4]: As Κ (0.432 N/m) is much larger than  0 (7 µN/m), it is practically determined by the   +   contribution in the model.For simplicity in the model,   and   were set to equal values.

Volume constraint
The volume incompressibility is enforced to ensure the RBC volume is conserved.This is achieved by an isotropic penalty term similar to the area incompressibility energy function: where Ω is the current cytoplasmic volume of the RBC, Ω 0 is the original cytoplasmic volume based on the mean corpuscular volume (MCV) of 100 fL and  Ω is the volume correction penalty coefficient set as 220 N/m 3 .

Section G: Plasma and cytosol transport model: Lattice Boltzmann Method (LBM)
The LBM solves the Boltzmann equation numerically by discretizing the mass and momentum of fluid particles on a fixed computational lattice.In the LBM, the behavior of the macro-system (defined by pressure, viscosity and density) is given by the kinetics and sum of its microstates.Microstates are defined in the discrete packets by the density distribution function (  ()) and the discretized transport equation involves the streaming and collision of   () in the lattice.In the streaming process, 19 discrete streaming directions and lattice velocities are defined as follows [7]: where ∆ is the space step size in the lattice, ∆ is the time step size, and symbol i denotes the flux direction index.The LB equation with a general force term can be expressed as [8]: where   is the body force term represented on the microstate and  is a relaxation parameter towards the equilibrium distribution (   ( , )) which can be obtained through: where  is the fluid density and  is the velocity; the two are calculated by the equivalency of the macroscopic state to the sum of microstates:  = ∑    ,  = ∑      ⁄  , and   = ∆ ∆ √3 ⁄ (speed of sound in the lattice).The weights,   , used in Eq (G3) are given by  0 = The microstate body force term in Eq (G2) is given by: where   is the fluid-structure interaction body force acting on the fluid from the IBM coupling with the RBC membrane solver.Proof of the equivalency between the LBM discretization and the continuum Navier-Stokes equations has been given through the Chapman-Enskog expansion of the LB equation [9] The pressure () and kinematic viscosity () of the fluid is represented in the LBM system by:  =    2 ;  =

Section C :
Fig B. RBC velocity data in the caudal artery CA (i), caudal vein CV (ii) and intersegmental vessels ISVs (iii) in trunk region pertaining to the simulation model domain (posterior to the yolk extension) for WT zebrafish at 2dpf and their corresponding RBC counts in the ISVs (iv).Orange grid lines in velocity plots (i to iii) indicate the time positions of systolic peaks in the CA velocity and show the phase lag between CA peaks and CV peaks.

Fig C .
Fig C. RBC velocity data in the caudal artery CA (i), caudal vein CV (ii) and intersegmental vessels ISVs (iii) in trunk region pertaining to the simulation model domain (posterior to the yolk extension) for Marcksl1a and Marcks1b double knock-out (referred to as Marcksl1KO in main text) zebrafish exhibiting indistinguishable flow and morphological difference from WT at 2dpf.Their corresponding RBC counts in the ISVs are shown in (iv).Orange grid lines in velocity plots (i to iii) indicate the time positions of systolic peaks in the CA velocity and show the phase lag between CA peaks and CV peaks.

Fig D .
Fig D. RBC velocity data in the caudal artery CA (i), caudal vein CV (ii) and intersegmental vessels ISVs (iii) in trunk region pertaining to the simulation model domain (posterior to the yolk extension) for Marcksl1a and Marcks1b double knock-out (referred to as Marcksl1KO in main text) zebrafish exhibiting moderate lumen diameter and RBC perfusion reduction phenotyping at 2dpf.Their corresponding RBC counts in the ISVs are shown in (iv).Orange grid lines in velocity plots (i to iii) indicate the time positions of systolic peaks in the CA velocity and show the phase lag between CA peaks and CV peaks.

Fig E .
Fig E. RBC velocity data in the caudal artery CA (i), caudal vein CV (ii) and intersegmental vessels ISVs (iii) in trunk region pertaining to the simulation model domain (posterior to the yolk extension) for Marcksl1a and Marcks1b double knock-out (referred to as Marcksl1KO in main text) zebrafish exhibiting extensive lumen diameter and RBC perfusion reduction phenotyping at 2dpf.Their corresponding RBC counts in the ISVs are shown in (iv).Orange grid lines in velocity plots (i to iii) indicate the time positions of systolic peaks in the CA velocity and show the phase lag between CA peaks and CV peaks.

Section D :
Fig F. RBC flow rate and blood flow rate comparisons between smooth geometry models (SGMs) for the caudal artery (CA) and caudal vein (CV) vessels.The results indicate the same flow-rate condition setting in the CA and CV amongst the 3 SGMs: 2.13, 2.14, 2.15 nL/min for RBC flow rate and 11.11, 11.04, 11.33 nL/min for blood flow rate in the CA for SGM1, SGM2 and SGM3; 2.07, 2.12, 2.16 nL/min for RBC flow rate and 12.45, 12.48, 12.44 nL/min for blood flow rate in the CV for SGM1, SGM2 and SGM3,.

Fig G .
Fig G. Pressure distributions for the 3 SGMs (i) indicating pressure gradient accentuation in ventral regions of aISVs in SGM1 and SGM2 due to the narrowing of lumen cross-sections in the dorsal to ventral direction , in contrast to smooth pressure gradients in aISVs of SGM3 that has constant lumen diameters throughout the vessel length (ii) ( +   ∆,  + ∆) −   (, ) = − |  (, ) −    (, ) Fig I.The velocity prediction of the LBM solver for straight circular tube of 30 µm and 150 µm in length, performed at different grid resolutions and different relaxation parameter (τ) values.(i) Velocity map fills of the lumen cross-section and (ii and iii) parabolic velocity profiles in the dorsal to ventral direction of the cross-sectional lumen plane at different LBM grid and τ settings.τ settings in (ii) achievable by density scaling in the LBM leads to better flow prediction accuracy as LBM grid is improved from coarse level to very fine level.When τ was allowed to change by not performing density scaling in (iii), the increase in grid resolution had contrary effects to solution accuracy.This highlights τ as an important parameter controlling the level of numerical diffusion at boundary condition regions seen in (iii) compared to (ii).In our main paper the τ value was fixed at 1.22.

Fig J.
Fig J.Simulation model sensitivity to the Dirac-delta stencil width () employed in the Immersed Boundary Method (IBM) fluid structure interaction solver.(i) RBC velocity maps for the wild-type (WT) network model performed with  = 1 and (ii) the resulting collision instability scenario where an RBC is artificially stuck on a lumen wall structure due to zero velocity on the near-wall grid.(iii) RBC velocity maps for the wild-type (WT) network model performed with  = 2 and (iv) the resulting robust collision resolved scenario where an RBC creeps away from the lumen wall protrusion thereby avoiding extensive deformation: as the stencil is wide enough to interpolate finite velocity to the collided RBC mesh nodes, the RBC will not be subjected to prolonged stretching that can break the code.(v) Velocity distribution histograms for  = 1 and  = 2 models for the RBCs flowing in the caudal artery (CA).(vi) The mean RBC velocity in the CA for  = 1 and  =2 models, indicating a 12% higher effective viscosity in  = 2 model as a result of a more diffusive stencil.