Angiogenic Microvascular Wall Shear Stress Patterns Revealed Through Three-dimensional Red Blood Cell Resolved Modeling

Abstract The wall shear stress (WSS) exerted by blood flowing through microvascular capillaries is an established driver of new blood vessel growth, or angiogenesis. Such adaptations are central to many physiological processes in both health and disease, yet three-dimensional (3D) WSS characteristics in real angiogenic microvascular networks are largely unknown. This marks a major knowledge gap because angiogenesis, naturally, is a 3D process. To advance current understanding, we model 3D red blood cells (RBCs) flowing through rat angiogenic microvascular networks using state-of-the-art simulation. The high-resolution fluid dynamics reveal 3D WSS patterns occurring at sub-endothelial cell (EC) scales that derive from distinct angiogenic morphologies, including microvascular loops and vessel tortuosity. We identify the existence of WSS hot and cold spots caused by angiogenic surface shapes and RBCs, and notably enhancement of low WSS regions by RBCs. Spatiotemporal characteristics further reveal how fluctuations follow timescales of RBC “footprints.” Altogether, this work provides a new conceptual framework for understanding how shear stress might regulate EC dynamics in vivo.


Introduction
Ang iog enesis, the process of growing new blood vessels from existing v essels, generall y occurs in r esponse to local stimuli.4][5] Three-dimensional (3D) WSS c har acteristics experienced due to red blood cells (RBCs) flowing through real ang iog enic micr ov ascular networks with complex surface topologies, howev er, ar e generall y unknown.Making such measurements in experiments is difficult, and in silico modeling that r esolv es the necessar y 3D geometric details and explicit RBC deformation is extr emel y challenging.Studies have shown intricate WSS patterns can arise due to complexities such as EC curv atur e 6 and 3D shapes, 7 and more recently due to complex 3D vessel morphologies and influence of RBCs. 8Such findings underscore the importance of elucidating 3D WSS complexities that occur during ang iog enesis and the associated micr ov ascular network remodeling in vivo, knowledge of which can help resear c hers build a new mec hanistic understanding of the processes in 3D leading to profound impacts affecting public health.
Ang iog enic processes are ubiquitous at all stages of animal life, and in both health and disease. 9The development and growth of organs in embryos is facilitated by ang iog enesis, 10 and continues up until adulthood. 9Processes associated with wound healing, [11][12][13] the menstrual cycle, 14 , 15 and pregnancy 16 , 17 inv olv e ang iog enesis, as do exercise 18 and muscle adaptation. 19 , 20The development and progression of many diseases also r el y on ang iog enesis. 21In cancer, the process enables invasion of tumors into tissue and can even supply a route for metastatic cells to exit and circulate to other regions of the body. 98][29] Other conditions such as psoriasis, 30 , 31 arthritis, 32 , 33 osteoporosis, 34 , 35 and macular de gener ation 36 , 37 all inv olv e ang iog enesis as well.
Current knowledge of hemodynamic c har acteristics that exist during ang iog enesis in vi v o has in large part been provided by in silico fluid dynamics modeling.For example, WSS v ariations ar ound isolated micr ov essel loop structur es in avian c hic k embry os w ere shown using a 2D single-phase blood model, and connections between specific WSS patterns and sprout locations were identified. 38Other work in retinal ang iog enic networks using a 3D single-phase flow model demonstrated general WSS values, which can lead to vessel regression. 39][47][48][49][50][51][52][53][54][55] While much has been learned regarding the WSS behavior associated with ang iog enic events, 3D aspects and influences of RBCs are still unclear.
Direct evidence that ECs respond to applied forces has mainly come from experiments focusing on isolated ECs.Specific WSS values eliciting EC responses in vitro , commonly measured in parallel plate flow-type devices, have been reported to range anywher e fr om 0.1 to 100 dyne/cm. 21 , 5 , 56-60While such findings importantl y esta b lish physiolo gical conte xt for WSS magnitudes and ang iog enic behavior, the significant rang e of WSS values spanning three orders of magnitude leads to questions concerning what is actually experienced by ECs during ang iog enesis in vi v o.Such questions are further complicated by findings from other in vitro works that have demonstrated strong EC responses to WSS spatial patterns or gradients, [61][62][63][64] as opposed to just single WSS values.A few recent in vivo works have demonstrated WSS values up to 2 dyne/cm 2 corresponding to EC rearrangement in developing retinal networks, 65 or during intussuscepti v e ang iog enesis in isc hemic muscle . 66Yet, what do 3D angiogenic WSS patterns look like in vi v o?This seemingly simple question r e pr esents a major gap in understanding because new v essel gr owth and adaptation occurs in vi v o thr ough 3D changes to complex surface morphologies, not as a 1D or 2D process or in response to a single WSS value.Moreover, since microcirculator y v essel diameters ar e similar to or smaller than individual RBC size, 67 , 68 RBC deformation and interactions coupled with vessel surface complexity strongly influence the fluid dynamics and near-wall behavior, as has been shown by recent highresolution modeling. 8 , 69This further emphasizes the importance of r ev ealing 3D ang iog enic WSS patterns because ang iogenic micr ov ascular networks hav e distinct geometric features such as looping v essel structur es, 38 fr equent br anc hing patterns, 70 and increased vessel tortuosity. 71o war d this end, we reveal and quantify in this work 3D WSS c har acteristics that occur in real ang iog enic microvascular networks by integrating high-fidelity RBC-r esolv ed sim ulations with imaging of the rat mesentery undergoing ang iog enesis.Specifically, we study the effects of ang iog enic remodeling on local WSS distributions that in-turn influence EC responses.Behavior is identified by first focusing on variations from vessel to vessel to establish context for describing highly localized 3D variations in both space and time in subsequent sections.Detailed maps of WSS spatial patterns are provided, which illustr ate ric h 3D beha vior influenced b y a wide range of ang iog enic geometrical features and RBCs.Altogether, the findings present a new 3D picture of the hemodynamic micr oenvir onment likel y experienced by ECs during ang iog enic network remodeling in vi v o and suggest a conceptual foundation for considering the r egulator y inv olv ement of shear str ess on EC biology.

Materials and Methods
The approach employed in the current work inv olv es integrating high-r esolution ima g es of ang iog enic micr ov ascular networks in the rat mesentery with high-fidelity RBC-resolved simulations in 3D.Output data from the simulations includes detailed information on the near-wall fluid dynamics, which were used to calculate 3D WSS and its spatial patterns.The sections below describe the methods employed for each of the r especti v e components.

RBC-Resolv ed Sim ulations
Three-dimensional computational fluid dynamics simulations are performed using an immersed boundary method-based approach for modeling biophysical flows in complex geometries. 72Salient details of the method are pro vided belo w, with a graphical overview provided in Figure 1 .The method has been extensi v el y v alidated a gainst both theor y and alternate n umerical methods, 72 as well as data from in vivo experiments. 8 , 69 , 73hat this n umerical appr oach can r e pr oduce what has been generall y observ ed in vi v o w as demonstr ated in prior w ork.Here , w e F igure 1. Overvie w of in silico method for modeling 3D RBC-r esolv ed flows in complex geometries.extend this to understand ang iog enic WSS c har acteristics based on real image data.
Blood is modeled as a suspension of 3D RBCs flowing with plasma, and the cells deform as they flow through the complex 3D micr ov essels.Each RBC is r e pr esented by an infinitesimally thin elastic membrane enclosing a viscous fluid r e pr esenting hemoglobin.Both plasma and hemoglobin are modeled as Newtonian fluids, with different dynamic viscosities; plasma pr operties ar e similar to w ater (1.2 cP), while the hemoglobin is modeled as being five times more viscous. 68Each RBC membrane is a surface defined by a Lagrangian mesh of 5120 triangular elements, and simulations each resolve approximately 1500 RBCs.The RBC resting shape is taken to be a biconcave discocyte with a major axis diameter of 7.8 μm, a surface area of 134.1 μm 2 , and a volume of 94.1 μm 3 , 68 , 74 which correspond to human RBCs.Differences in major diameter between human and rat RBCs are approximately 0.8 μm, 75 Note that this considers a dynamic viscosity field that varies in time and space, or μ = μ( x , t ) .At a gi v en time and instantaneous RBC configuration, the spatial locations outside RBCs have plasma viscosity, while inside RBCs have a viscosity five times higher.Mathematically we model this as μ( x , t ) = μ p + ( μ c − μ p ) I ( x , t ) , where μ p and μ c are viscosities as denoted in Figure 1 , and I ( x , t ) is an indicator function used to tr ac k the separate fluids.Full details on the model and implementation can be found in prior works. 72 , 73The governing fluid flow equations are solved for the fluid velocity ( u ) and pressure ( P ) on a fixed, uniform Eulerian mesh with a projection method using a finite volume/spectral appr oach, wher e the body for ce term F incorpor ates stresses fr om the deforma b le cell model described below.Grid spacing for the fluid field is appr oximatel y 200 nm.
Complex vascular walls are modeled with a sharp-interface ghost node (GN) IBM, which decomposes the Eulerian computational domain into a fluid domain inside the vessels and a solid domain outside.The interface conditions at the walls are enforced when solving the fluid flow equations via constraints imposed at the Eulerian mesh points immediately outside the fluid domain (points termed as ghost nodes).Each GN has a corr esponding ima ge point (IP) in the fluid domain, which is the mirr or ima ge of the GN acr oss the near est location on the wall boundar y, termed boundar y interce pt (BI).The velocity at the BI point is taken as the av era ge of the GN and IP velocities, or . Thus, to enforce a no-slip condition at the BI point, the velocity at the GN is written as u GN = −u IP for u BI = 0 .While the GN velocity is necessarily at an Eulerian mesh point, the IP is not.Thus, we r e pr esent the velocity at the IP based on the velocities at the 8 mesh points surrounding the IP using a trilinear interpolant of the form u IP = 8 i= 1 β i u i , where β i are the interpolation weights.
The membrane stresses in response to deformation are calculated using the finite element method, with loop elements used as a subdivision surface for the force calculations. 76 , 77The surface force density f e at each triangular element is determined b y e valuating the surface gradient of the Cauchy tension tensor ( T mn ), and membrane ( σ mn ) and bending ( ν mn ) stresses as f e = ∇ s ( T mn + σ mn + ν mn ) .Force densities are transferred from elements to vertices using box-spline shape functions, and the vertex forces f v are spread to the Eulerian mesh via F using a delta function in terms of Eulerian grid location x and membrane vertex location X : This same delta function facilitates interpolation of the Eulerian velocity field to update the membrane vertex locations and deform the cell.To evaluate the Cauchy tension tensor an energy function developed for a RBC 78 is used to determine the str ess r esponse to the strain as , where G s is the membrane shear elasticity, and I 1 , I 2 are the strain invariants.The membrane can undergo significant deformation yet is nearly ar ea-incompr essib le, in a gr eement with physiological behavior of RBCs. 79Resistance to surface area changes is enforced by the constant C in the previous equation for strain energy.A larger value of C results in increased area preservation but at the expense of numerical stability.Here , w e have used C = 100, whic h w e observ e kee ps the surface ar ea changes to less than 0.1%.The Helfrich bending energy is used for the membrane and bending stresses, written as σ mn = k b 2 ( 4 κ 2 A mn − 8 κ B mn ) and ν mn = k b 2 ( 4 κ A mn ) , r especti v el y, wher e k b is the bending modulus, κ is the mean curv atur e, and A and B are the metric and curv atur e tensors.Computations follow that described elsewhere. 80imulation snapshots shown in Figure 2 depict the entire domain of each modeled network, with vessels at the borders being either inlets or outlets.Flow rates are imposed at the inlets and outlets to achieve physiological effecti v e shear rates generally in line with reported values in the literature for the mesenteric microcirculation. 81For each network, three different flow conditions are modeled, and denoted as baseline, 0.5x and 2x.Baseline shear rates at boundaries range from approximately 30 to 500 s −1 , and are scaled for the other flow cases based on the denoted factor.Details are provided in the Supplementary Materials.At inlet vessels, hematocrit levels are maintained by adding RBCs from external, linked simulations involving RBCs flowing through straight tubes of the same diameter as the r especti v e inlets.The inlet hematocrit values are prescribed generally based on vessel diameter following physiological values reported for the microcirculation, 81 , 82 and range from appr oximatel y 12% to 36%.Details are provided in the Supplementary Materials on how RBC injection is handled at inlets, as well as the initial placement of RBCs throughout the networks.It is noted that for some inlet vessels the hematocrit slowly decr eases ov er time, and this behavior is discussed in the Supplementary Materials.

Calculation of WSS and WSSG
The high-resolution 3D velocity field enables calculation of the WSS and its gradient (WSSG) on the complex micr ov ascular w all surfaces.To facilitate study of these quantities, we calculate the traction vector ( t ) at each vertex of the mesh defining the vascular walls.Specifically, t gives the force per area at the surface, and is related to the stress tensor σ by t = n • σ , where n is the unit normal vector at the surface.Following previous work, 8 we use a local cylindrical coordinate system for the WSS calculation at each wall point due to the geometric complexity.The axial direction ( e a ) points along the direction of the mean undisturbed flow near the wall, based on the idea that the mean flow moves in the axial direction.The radial direction ( e r ) is in the inw ard-normal dir ection, and the circumferential component ( e θ ) is defined as e θ = e a × e r .This approach provides a consistent means of calculating WSS at all wall points and, furthermor e, pr ovides a conv ection that is similar in spirit to that used with standard Poiseuille flows in straight tubes.With this convention, the viscous components of the traction vector can be r e pr esented as t a = μ p ∂u a ∂r , t θ = μ p ∂u θ ∂r , and t r = 2 μ p ∂u r ∂r , where μ p is the plasma viscosity since the calculation is at the vessel wall.While we observe the axial component to be many orders of ma gnitude gr eater than the other components, the WSS values r e ported her e ar e taken to be the magnitude of the traction vector to include all data.
To calculate the WSS, the deri v ates in the equations for the traction components ar e n umericall y ev aluated using secondorder differ encing inw ards fr om the surface.The v elocity components at each location in the stencil are first interpolated from the velocity field data, and then converted to the local coordinate system at the vertex.Thus, given a 3D velocity field determined from an RBC-resolved simulation, the resulting WSS at all surface points is calculated following this pr ocedur e. Subsequently, with the WSS determined its spatial gradients can be calculated.We calculate the WSSG components in both the axial and circumfer ential dir ections thr ough surface interpolation of the WSS field and using second-order central differencing.

Imaging of Angiogenic Micro v ascular Netw orks and Digital Reconstruction
Netw ork images w ere obtained by re-imaging of angiogenic adult rat mesenteric tissues from a previous study. 83As described therein, ang iog enesis w as stim ulated via i.p. injections of compound 48/80, a mast cell de gr anulator.4][85][86] The representative ang iog enic network images for the current study were obtained from adult Wistar rat mesenteric tissues 10 d post stimulation. 86Our previous c har acterization of micr ov ascular r emodeling r esponses post stimulation in these tissues documented that vascular area per tissue is increased by day 10.Capillary sprouts per vascular ar ea wer e incr eased by 3 d post the injection protocol and r emained incr eased at day 10.At day 10, vascular length per vascular ar ea w as incr eased.Remodeled micr ov ascular networks also displayed an increase in v en ule tortuosity and capillary plexus regions reflective of patterning changes.
Post stimulation, perfusion of fixa b le 40 kDa FITC-dextran identified vessel lumens.Vessel segments wer e additionall y confirmed by labeling of blood ECs for PECAM.Mesenteric tissues were whole mounted and observed en face.The mesentery's thinness (20-40 μm) ena b les ima ging of intact networks at the vessel and cellular levels and eliminates the need for tissue sectioning.Micr ov ascular networks exhibited increased segment density, 83 which is characteristic of ang iog enesis regardless of the stimulation method.Thus the 48/80 stimulated networks are sufficient to provide examples of ang iog enic network re gions.Tw o representative netw orks of high-density capillary re gions w ere selected and captured using 10x (dry, NA = 0.3) objecti v es on an inv erted micr oscope (Nikon Ti2) coupled with an ANDOR Zyla sCMOS camera.
For the present work, the ang iog enic microvascular regions ( F igure 2 ) w ere digitally reconstructed for 3D RBC-resolved simulations and anal ysis.PtcCr eo 87 w as used to build the networks by first importing images on a plane, a reasonable approximation gi v en the mesenter y's thinness as noted a bov e , and extr acting approximations for vessel centerline trajectories.Any outof-plane vessels are identified through visual inspection, and centerlines modified accordingly.The 3D interconnected vessel structur es ar e then cr eated by swee ping circular cr oss-sections along the trajectories, with sizes adjusted locally in accordance with the images.The shapes of connection regions are similarly patterned to match the images.It is noted that the use of centerlines and locally circular lumen shapes is a modeling assumption due to images being 2D.Geometric details of the 3D network surfaces are then exported in .stpformat and imported into the Gmsh softw ar e 88 for surface triangulation.Such meshing r esults in highly resolved 3D discretized microvascular surfaces that comprise eac h netw ork, with triangles having edge lengths of appr oximatel y 150 nm.The ang iog enic network in Figure 2 A-C has 9.6 million vertices and 19.2 million triangular faces, while that in Figure 2 D-F has 7 million vertices and 14 million triangular faces.Geometries are stored in .stlformat, which are used to r e pr esent the micr ov ascular w alls for the RBC-r esolv ed sim ulations.

Angiogenic WSS Patterns Exhibit Significant Heter ogeneity fr om V essel to V essel
Surface contours of the time av era ged WSS (TAWSS) pr edicted for each of the ang iog enic micr ov ascular networks ar e shown in Figure 2 from representative simulations.In addition to the TAWSS maps, Figure 2 A-C and D-F provide snapshots of RBCs flowing through each respective network along with the image data from which the networks were reconstructed.Movies are also provided in the Supplementary Materials.Time-averages are performed over the entire simulation time for all cases, which is on the order of one second.
Significant TAWSS heter ogeneity acr oss eac h netw ork is evident from the representative data shown in Figure 2 C, F, where WSS contour values can be qualitatively observed to v ar y by almost two orders of magnitude between different vessels (see Figure 3 for zoomed-in perspectives of representative regions).To quantify this observation, the distribution of predicted TAWSS on a per-vessel basis (ie, spatiall y av era ged ov er all points in a vessel) is gi v en in Figur e 2 G for all simulations as a function of vessel diameter.From this figure, we observe the magnitude of the vessel-to-vessel TAWSS variability to be generall y de pendent on v essel diameter.Variations ar e most significant in the smallest capillar y v essels (3-15 μm diameter), wher e v alues range fr om 0.01 d yne/cm 2 to 130 d yne/cm 2 .With increasing diameter the vessel-to-vessel variability decreases, with TAWSS values approximately ranging from 0.01 dyne/cm 2 to 40 dyne/cm 2 in vessels greater than 25 μm in diameter.Additional perspecti v e on heter ogeneity on a per-v essel basis is pr ovided by the heatmap in Figure 2 H, which gi v es the TAWSS per vessel in terms of both vessel length and diameter.This data show that both high and low TAWSS vessels are scattered across all vessel lengths for diameters less than 20 μm.Similarly, for larger diameters the TAWSS per-vessel values (albeit lower in magnitude) are scattered across all vessel lengths without a clear pattern.These findings suggest that the TAWSS varia bility fr om v essel to v essel is not str ongl y influenced by v essel length in the reconstructed ang iog enic networks, only the diameter ( Figure 2 G).
We observe that the enhanced TAWSS heterogeneity in small diameter vessels to be str ongl y influenced by the presence of RBCs.This is caused by increased confinement leading to increased RBC deformation and wall interactions, which, perhaps expectedly, leads to the high TAWSS values in small diameter vessels.This beha vior, howe ver, also contributes to low TAWSS in such vessels.The simulations predict per-vessel TAWSS values to occur in direct proportion to predicted vessel shear rates (see Supplementary Material), and in a few vessels ther e ar e negligib le flows, whic h natur ally leads to ne gligible TAWSS values.Some of these lo w-flo w vessels have diameters 3 μm or less, which through physical occlusion effecti v el y pr event the passage of RBCs.Through netw ork structure , cells are routed to other adjacent vessels, which in many cases incur relati v el y high TAWSS.Inter estingl y, we also observe small diameter v essels dev oid of RBCs can hav e TAWSS v alues upw ards of 40 dyne/cm 2 (see Supplementary Material and Discussion).Collecti v el y, these behaviors are only captured by explicitly resolving the RBCs.We observe these patterns of vessels with low TAWSS adjacent to vessels with high TAWSS to occur for either single vessels or groups of vessels.Representative examples are provided in Figure 3 .

Complex 3D Spatial Variations Occur at Sub-EC Length Scales
Close inspection of the TAWSS surface contours in Figures 2 and 3 r ev eals complex 3D spatial patterns underlying the perv essel TAWSS v alues.Re pr esentati v e examples of these patterns ar e pr ovided in Figur e 4 for one of the ang iog enic network regions.
As evident from the representative behavior depicted in Figure 4 A, TAWSS variations are observed on surfaces within indi vidual v essels.Namel y, both the ma gnitudes of the localized TAWSS values experienced at surface points as well as the degree of spatial variation among the vessels shown are highl y heter ogeneous.The 3D surface contours in Figure 4 B provide a zoomed-in perspecti v e illustr ating c har acteristic behavior we observe, namely the formation of TAWSS hot and cold spots within vessels and sharp changes over relatively small distances.For the example vessel in Figure 4 B, based on the TAWSS spatial patterns and the scale bar shown, significant TAWSS variations occur over distances that are O ( μm), which are shorter distances than typical size of ECs that line vessel walls.
To quantify the magnitude of localized 3D spatial variations within vessels, we first compute the standard deviation ( σ ) of the TAWSS among all surface points in each vessel.Data shown in Figure 4 E provides the predicted distribution of σ versus vessel diameter for all vessels.In the smallest capillaries σ values can reach as high as 60 dyne/cm 2 , while negligible values ar e observ ed as well in a few vessels.This date indicates that within vessels of similar diameter the degree of spatial variation in TAWSS can range over an order of magnitude.Spatial variations within vessels of larger diameter among the data points in Figure 4 E show σ values reaching approximately 20 dyne/cm 2 , with a lower degree of vessel-to-vessel variation as compared to the smaller diameter vessels.
The structure to the TAWSS spatial patterns in Figure 4 B qualitati v el y exhibits a dir ectional pr efer ence, meaning sharp changes can occur in specific directions with respect to vessel orientation and flow direction.For example, the hot and cold spots associated with the identifying arrows in Figure 4 B occur in sequence along the mean flow direction (see Figure 4 A for flow dir ection), wher e a TAWSS hot spot is followed by a cold spot.This indicates that a dominant spatial change to the TAWSS occurs here in the direction of the mean flow or the axial direction.For this example, this occurs due to a temporary increase in vessel area based on the local morphology, which causes the observed spatial patterns exhibiting a negati v e TAWSS gradient in the axial direction (ie, a region of relatively high TAWSS transitioning to a region of low TAWSS).From the high-resolution simulation output we can mathematically quantify such behavior by computing the gradient of the TAWSS field in the axial direction (TAWSSG axial ).Contours of this quantity are shown in Figure 4 E calculated from the TAWSS field shown in Figure 4 B. Figure 4 D pr ovides TAWSSG axial av era ged ov er all points in each v essel v ersus vessel diameter.Here, the absolute values are plotted to generall y r eflect the ma gnitude of the pr edicted spatial v ariations.
To add further perspecti v e to the spatial variations quantified as spatial gradients, we also compute the TAWSS gradient in the circumfer ential dir ection (TAWSSG circ ).Contours of this quantity are shown in Figure 4 F, which augments the three-dimensional nature of the observed TAWSS spatial variations in the above example.Namely, the TAWSS does not simpl y decr ease uniformly as the vessel cross-sectional area locally increases, but r ather c hanges in a muc h more complex and directionally biased manner.Pronounced TAWSS axial gradients also occur within r egions wher e v essels connect to one another, as in the lowermost horizontal vessel in Figure 4 E. For this example, as the flow enters the smaller diameter vessel the fluid accelerates resulting in incr eased w all shear rates and incr eased TAWSS in the axial direction.As with the previous example , how ever, this spatial variation is much more complex than a simple increase in value.Instead, TAWSS values increase and spatial patterns c hange .This is quantified in Figure 4 E, F by the regions of pronounced TAWSSG axial and circumferential components.
Characteristic ang iog enic network structur es ar e micr ov ascular "loops," where three or more vessels are connected to one another in an arrangement r esemb ling a loop structure.Re pr esentati v e examples ar e pr ovided in Figure 5 A. From our observation of large TAWSS spatial changes at locations where vessels connect to one another, we also analyzed behavior in these loop structures.A hemodynamic behavior we observe that affects the TAWSS is at least one of the vessels comprising a loop typically has negligible flow of RBCs.For the ang iog enic networks considered in this w ork, w e identified 25 micr ov ascular loop structures and among them 17 exhibited this behavior (see Supplementary Material).Due to such flow discr e pancies amongst connected vessels, our simulations reveal that the resulting TAWSS patterns can hav e orders-of-ma gnitude spatial variations over v er y small distances.Examples are shown in Figure 5 A, where for these loop structures one of the connected br anc hes is a lowflow (or negligible flow br anc h).Inter estingl y, we also r e port low flow br anc hes in whic h the flow dir ection v aries ov er time, an example of which is gi v en in Figur e 5 A and denoted by the yello w arro ws.This occurs her e just downstr eam of a bifurcation, in a cr oss-v essel wher e we observ e RBCs to flow v er y slowl y and in either dir ection de pending on the local time-dependent configuration of flowing RBCs and their interactions (movie provided in the Supplementary Material).
Gi v en the significant degree of TAWSS spatial variations within vessels, we also analyzed the tendency of individual vessel surfaces to experience high or low TAWSS values.To achieve this, we considered TAWSS values across all surface points for a vessel, and first binned the TAWSS across the range of values experienced for that vessel.Bin values for each vessel span the full range for the v essel fr om low to high TAWSS and are relative to the minimum and maximum TAWSS values for each vessel.Values are distributed over nine bins for each vessel, so that the lowest and highest bins correspond to the low and high TAWSS values for the vessel.Then, for each bin, we determined the corr esponding v essel surface ar ea fraction to provide a histogram of TAWSS behavior on a per-area basis for each vessel.These data ar e pr ovided in Figur e 5 B, which gi v e the av era ge ov er all v essels.Inter estingl y, we observe that on aver age , the area fraction of a vessel experience high TAWSS tends to be v er y small, wher eas areas experiencing low TAWSS tend to be large.

RBCs Alter T A WSS Behavior in Angiogenic Networks
The presence of RBCs expectedly alters the TAWSS over that experienced based on pure plasma flow alone.We quantify this beha vior b y comparing the per-v essel TAWSS fr om the RBCr esolv ed sim ulations to the per-v essel TAWSS fr om sim ulations in the ang iog enic networks without any RBCs (ie, only plasma).Differences manifesting in various aspects of the TAWSS character are quantified in Figures 6 and 7 .
Comparati v e differ ences between TAWSS with and without RBCs are shown in Figure 6 A. For each vessel, we compute the mean, maximum, minimum, and standard deviation of the TAWSS among all surface points and plot the ratio of these values v ersus v essel diameter for sim ulations with and without RBCs.Quantities are separately plotted in Figures 6 A-i and Aii, in order to isolate differences in behavior we observe.That is, in Figure 6 A-i, the mean, maximum, and standard deviation ratios are plotted, and these quantities all exhibit similar behavior; values have similar magnitudes with respect to diameter, and are greater than unity for most vessels (ie, RBCs increase the quantity being considered).Specific values are gener ally betw een 1.5 and 3 for most vessels, while in some vessels v alues ar e less than unity.Reasons for the later ar e discussed in a subsequent section on local spatial variations, but in general occur due to vessels with enhanced cold spots caused by RBCs and/or vessels with very short length and complex morphology.The behavior in Figure 6 A-i contrasts that of the ratio for minimum TAWSS , whic h is shown in F igur e 6 A-ii.While many v alues are indeed greater than unity, interestingly many vessels also hav e v alues less than unity.To best illustrate this, values are plotted in Figure 6 A on a log scale.As indicated in the figure, data points that fall below the 1.0 line r e pr esent v essels in which the presence of RBCs has reduced the minimum value of TAWSS experienced in the vessel or enhanced the low TAWSS c har acter.
The influence of RBCs on the TAWSS will naturally lead to influences on the TAWSSG.To quantify this influence, we plot the ratio of TAWSSG with and without RBCs in Figure 6 B.Here , w e denote this r atio as α = TAWSS G RBC / TAWSS G plasma to facilitate comparing the roles of the directional components via Figure 6 B-i and B-ii.The α values for the axial and circumfer ential components ar e plotted in Figur e 6 B-i v ersus v essel diameter.As shown, the RBCs generally increase the TAWSSG acr oss all v essel diameters, will v alues incr easing by as m uch as five times compared to the plasma simulations.An interesting result is the increased α values for the circumferential component as compared to the axial component.That is, the presence of RBCs significantly enhances the TAWSSG circumferential component compared to the axial component.Recent work r e ported similar behavior, 8 which we also show here for the ang iog enic networks considered.To quantify the magnitude of enhancement to the TAWSSG circumferential component compared to the axial component, we plot the ratio of these α components in Figure 6 B-ii.This generally shows that the enhancement ranges from 2 to 5 times across the range of vessel diameters.
For the ang iog enic networks considered, simulations were performed over a range of driving flow rates (or, shear rates) at network boundaries, to inv estigate de pendence of r esults on imposed boundary conditions.This is especially due to the RBCr esolv ed flows her e, and the potential for a non-linear relation between the perfusion distribution across the networks and rate of flow increase at the boundaries due to the presence of the deforma b le RBCs.To inv estigate this r elation, we consider data from our baseline simulation (see Materials and Methods) and compar e a gainst data fr om sim ulations wher e flows at all  boundaries have been modified by factors of 0.5 and 2.0 (r eferr ed to here as γ = 0.5 and γ = 2.0).In Figure 7 A, we first show the influence on the distribution of TAWSS v alues acr oss the networks.As can be seen, the mean v alues r oughl y incr ease in pr oportion to γ .The TAWSS v aria bility acr oss the networks also incr eases with incr easing γ .To better understand the TAWSS values that cause this behavior, in Figure 7 B, we plot the ratio of TAWSS at each γ to the baseline case.Inter estingl y, for v essels with diameters less than 10 μm, the TAWSS per vessel does not scale in proportion to γ as would occur in the absence of RBCs (see right-most plot in Figure 7 B for behavior from the plasmaonl y sim ulations).Instead, we observ e highl y scatter ed behavior, which emphasizes the role played by RBCs in influencing the perfusion distribution and resulting TAWSS behaviors.

RBCs Enhance Both Hot and Cold Spots
TAWSS hot and cold spots, defined as localized and concentrated regions of high or low TAWSS values, are predicted in most vessels throughout the networks.The high and low values that define spots ar e r elati v e quantities that depend on local TAWSS magnitudes, whic h w e qualitati v el y identify by adjusting local contour values to the TAWSS range experienced in the field of view for network sections.While specific TAWSS values defining a hot or cold spot will necessarily vary, local variations by r oughl y 10 dyne/cm 2 / μm or greater tend to form spots that can be visually identified.
Hot or cold spot formation is necessarily connected to highly localized TAWSS spatial gradients, which conce ptuall y follows the previous discussion associated with Figure 4 B. The enhancement of TAWSS hot spots by RBCs is expected, and we find this behavior can manifest by increasing both the magnitude and area of the hot spot.For TAWSS cold spots, enhancement inter estingl y causes the minimum TAWSS to be even smaller and the region much more focused.Such behavior underlies the behavior quantified in Figure 6 B-ii.Representative examples of hot and cold spot patterns and influence by RBCs are provided in Figure 8 , which complement the patterns in Figures 4 and 5 .
A primar y location wher e RBCs enhance TAWSS hot spots is within regions where vessels connect to one another.Depending on the local geometry, however, hot spot enhancements can manifest differ entl y.For the r e gion shown in F igure 8 A, the RBCs enter from the upper left and wrap around the pillar-like obstruction as they flow downstream.For this case, the nearwall motion of the RBCs causes a TAWSS hot spot on the side of the pillar as noted by the red arrow in the figure , whic h is enhanced compared to behavior for plasma-only simulations by nearl y fiv e times (see Supplementar y Material).Such enhancement is intuiti v e gi v en the expected increase in near-wall velocity gradient caused by the presence of RBCs.In contrast, other behaviors shown in F igure 8 A demonstr ates hot spot formation in non-intuiti v e locations.This can be seen in the lower-right region of the figure, where due to increased confinement, an RBC will tend to linger on the bifurcation apex and temporarily obstruct the flow.Other RBCs and plasma are subsequently r e-r outed locall y ar ound the b locka ge and cause incr eased v elocity gradients through these temporarily reduced flow areas.This leads to TAWSS hot spots on the outside wall, as noted by the red arrows.The formation of these hot spots is due to the discrete, particulate nature of the 3D RBCs, and would not occur with a contin uum b lood model.Along these lines, we note that such behavior is not enhanced by RBCs as is observed for other types of hot spots, but rather is caused by RBCs.
This bifurcation apex in Figure 8 A marks the location of a TAWSS cold spot, as indicated by the b lue arr ow in the figure .Suc h a cold spot is also observed in the plasma-only simulations and occurs in both cases because it is a flow stagnation point.The impact of RBCs lingering on the apex interestingl y r educes the TAWSS ma gnitude compar ed to the plasmaonl y sim ulations, but also causes the area of the cold spot to be reduced.Such cold spot enhancement, which manifests in a mor e focused r eg ion, is g i v en in Figur e 8 B, wher e TAWSS contours are shown for both RBC and plasma-only simulations.The reduced cold spot area is generally caused by the narrowing of the cold spot due to adjacent hot spots enhanced by the near-w all mov ement of RBCs flowing thr ough the bifurcation and moving downstream.This adjacent hot spot enhancement is denoted in Figure 8 B by the red arrows in the bifurcation on either side of the cold spot, along with comparison to the contours from the plasma simulation.The existence of the hot spot immediately to the left of this cold spot in the bifurcation is caused by changes to the local surface curv atur e of the v essel w all, which is another primar y means by whic h w e observe hot and cold spots to form.At this location, the vessel bends back and forth in the do wnstream flo w direction, resulting in a region of high surface curvature at the location of the hot spot, followed by a region of low surface curv atur e wher e the TAWSS is reduced.Such variations in TAWSS following the v essel curv atur e generall y occur due to the shifted v elocity pr ofile for cr ee ping flows in curv ed v essels; the velocity profile in a curved tube at low Reynolds numbers will be skewed to war d the side of the vessel with the higher curvature resulting in larger near-w all v elocity gradients. 8 , 69For the pr esent work in angiogenic vessels, this manifests in patterns such as that shown in Figure 8 B within the vessel feeding the bifurcation where hot spots shift back and forth between sides following changes to the surface curv atur e.This behavior is also observed for the corresponding TAWSS contours for plasma simulations, demonstrating that this phenomenon is caused by the geometry but enhanced by the RBCs.This can also result in TAWSS cold spots such as that denoted by the blue arrow on the left-most side of Figure 8 A.

Time-Dependent WSS Fluctuations Follow Timescales of RBC "Footprints"
The presence of individual deformable RBCs flowing through small micr ov essels cr eates WSS spatial patterns that mov e with the RBCs.In the present work, this time-dependent behavior occurs entir el y because of the RBCs, as the boundary conditions are constant.This enables isolation of the RBC contribution to the time-dependent WSS behavior, which occurs on different timescales than the car diac c ycle.For example, an RBC flowing through a microvessel at a velocity of 0.5 mm/sec will traverse a distance of 10 μm in 0.02 s.From the simulation data, we observe the resulting RBC "footprint" on the WSS, as termed in previous work in idealized vessels, 89 to appr oximatel y follow such timescales resulting in WSS fluctuations.Figure 9 provides To illustrate the relationship between instantaneous WSS spatial patterns and the movement of individual RBCs, the inset to Figure 9 A depicts WSS contours at three sequential time points along with the corresponding RBC positions.At each time point, the WSS contours exhibit spatial variations, and these can be observed to change following the RBC motion.This observation augments the concept of WSS heterogeneity by revealing how the sub-EC scale spatial variations are constantly in flux.To quantify the degree of WSS variations in both space and time , F igur e 9 B-D pr ovide example time signals for different WSS spatial statistics from two separate vessels identified in Figure 9 A. To re veal beha vior resulting from different hemodynamics expected in ang iog enic networks, we focus on (1) a small capillar y wher e RBCs flo w in single file and (2) a lar ger diameter v essel wher e RBCs flow in m ultifile fashion.The figur es focus on behavior over a time interval of 0.07 s in order to highlight the timescale of the fluctuations in quantities.
First, shown in Figure 9 B is the time signal for the mean WSS (ie, WSS spatiall y av era ged ov er all points in the vessel at each time instant).As can be seen, the small capillary fluctuates r oughl y between 45 and 55 dyne/cm 2 , with an oscillation timescale of appr oximatel y 0.021 s.We calculate this timescale as the inverse of the dominant frequency resulting from a fast Fourier transform of the time-series data.Temporal variations for the larger vessel are smaller, with values ranging fr om appr oximatel y 25 to 30 dyne/cm 2 .Next, in Figure 9 C, we plot the maximum vessel WSS versus time where the small capillary fluctuates between 120 and 190 dyne/cm 2 with an oscillation timescale of appr oximatel y 0.013 s, and the larger vessel from 60 to 130 dyne/cm 2 with a similar timescale.The differences between the magnitude of these values and their fluctuations are significantly enhanced compared to the mean fluctuations, providing a measure quantifying the degree to which the spatial v ariations ov er the v essels c hange with time .To quantify the time-evolving spatial variations further and more directly, we plot the WSS standard deviation over the two vessels versus time in Figure 9 D. The time signals for the standard deviation show values ranging from 15 to 26 dyne/cm 2 for the small capillary with an oscillation timescale of approximately 0.016 s, and 7 to 12 dyne/cm 2 for the larger vessel.To fully quantify oscillation timescales for the WSS metrics considered (ie, mean, maxim um, minim um, and standard deviation), data for all vessels is provided in the Supplementary Material.Similar to other quantities in previous sections, the oscillation timescales show significant heterogeneity in smaller capillaries, which decreases with increasing diameter.
Having esta b lished and quantified WSS fluctuations in r e pr esentati v e ang iog enic v essels, we also quantified time-de pendent fluctuations pr edicted acr oss all network sim ulations.To facilitate this, we calculated for each vessel the root mean squared (RMS) of the fluctuations for WSS mean, maximum, minimum, and standard deviation time signals.The data are provided in Figure 10 as a function of vessel diameter.As with our previous observations on TAWSS, the degree of heterogeneity can be seen to generally increase with decreasing diameter.For the WSS mean and standard deviation trends, RMS values are very similar, and generally range from 0.5 to 30 dyne/cm 2 .In contrast, RMS values for the maximum WSS values in each vessel range from 10 to 500 dyne/cm 2 , indicating significant temporal fluctuations can occur in some capillaries.A similar range for RMS of the minimum WSS in each vessel is observed, except magnitudes from 0.01 to 10 dyne/cm 2 .Similarity between the mean and standard deviation RMS values occurs because the WSS spatial distributions tend to move with the RBCs as mentioned, and these quantities r e pr esent a ggr egate behavior ov er whole v essels.This results in more gradual change with time r elati v e to the RMS values of the WSS extremum, which are naturally much mor e v olatile.

Discussion
The findings in this study r ev eal high-r esolution 3D WSS patterns in ang iog enic micr ov ascular networks based on r eal image data.Specifically, the results provide novel quantification of 3D WSS magnitudes, spatial, and temporal characteristics based on 3D RBC-r esolv ed sim ulations in r eal ang iog enic network structures.While 3D WSS variations over sub-EC length scales have been previously suggested for micr ov ascular networks by prior works, for example, 7 , 8 , 90 these previous studies inv olv ed either idealized or unstimulated networks.The qualitati v e de pictions and quantitati v e anal yses pr esented her e thus provide insights to the potential shear stress microenvironment likely experienced by ECs in real ang iog enic scenarios.Consequently, the findings lay the groundwork for furthering mechanistic understanding of EC behavior in remodeled networks and emphasize the need to consider shear stress as non-uniform stim uli that v aries ov er space and time in a highly localized manner.
An underlying goal of this work is to provide representati v e insights and quantifications of the types of WSS patterns that can occur when 3D RBCs and ang iog enic network morphologies are resolved.With this resolution, and the significant degree of WSS spatial variability shown here to result, questions then arise as to how such WSS complexity could possib l y guide ang iog enesis.While the current understanding is based on more coarse-grained representation, the current findings pr ov oke r e-consideration of the r oles pla yed b y WSS in affecting ang iog enic behavior.As models become mor e adv anced, the incr eased gran ularity pr ovides an opportunity to build new mechanistic models to understand ang iog enesis in a more localized and idiosyncratic fashion, even on a patient-specific basis.
The prediction of vessel-to-vessel T AWSS heterogeneity , and the significant range of values among different vessels of similar caliber, provides baseline behavior to esta b lish context for subsequent analysis.We have shown that on a per-vessel basis, TAWSS values can range anywhere from negligible to 130 dyne/cm 2 .In vitro works investigating EC responses to applied forces have shown values in this range to have ang iog enic significance.For example, works by dela Paz et al. 3 and Gloe et al. 4 demonstrated that when ECs were exposed to shear stresses ranging from 10 to 30 dyne/cm 2 the secretion of pro-ang iog enic growth factors (ie, bFGF) and the upregulation of growth factor r ece ptors (ie , VEGF-R2) w er e pr omoted.Similar str ess v alues were shown to promote EC sprouting, 1 while others have shown EC responses triggered by lower stress values, for example, 3-10 dyne/cm 2 , 58 , 91 or higher stress values upwards of 100 dyne/cm 2 . 59While such works have provided important data for inferring the force magnitudes experienced by ECs during angiogenesis, in vitro conditions are known to be very different from in vi v o envir onments. 92The dir ect calculations of TAWSS ma gnitudes r e ported her e per v essel fr om r eal ang iog enic networks ar e generall y in a gr eement with what can be inferr ed based on the in vitro works.However, our findings importantly show and quantify how r e pr esenting a vessel with one WSS value neglects the rich 3D spatial variations that are likely observed during ang iog enesis in vi v o.
Underlying the TAWSS vessel-to-vessel variability is heterogeneity in blood flow distributions throughout the networks.Nota b l y, we observ e behavior similar to recent work that r e ported heter ogeneity of RBC perfusion among connected v essels during micr ov ascular r emodeling. 93This alludes to the concept of a preferential path taken by RBCs through networks, resulting in regions where some vessels are always filled with RBCs and other vessels are always empty.The impact of this behavior on TAWSS naturally contributes to the TAWSS perv essel heter o geneity, and e xpectedly in our simulations we observ e v essels with negligib le TAWSS v alues to be generall y devoid of RBCs.Interestingly, however, the opposite is not always true; v essels generall y dev oid of RBCs do not al w a ys ha ve negligib le TAWSS v alues.In v essels with diameters generall y less than 5 μm, we observe that RBCs can, in some cases, be routed to other connected vessels due to mor e fav ora b le hydr odynamic resistance.Yet, the RBC-free vessel is non-negligibly perfused b y plasma.F or suc h cases, w e observe that TAWSS values can r each upw ar ds of 40 d yne/cm 2 (see Supplementary Material), which per the a bov e discussion is within a physiologically relevant r ange .This brings to light a potential limitation of inferring v essel WSS v alues based solel y on visual observ ation of flowing RBCs, particularly for smaller capillaries.
While single per-vessel TAWSS values can provide a first approximation of the EC microenvironment during ang iog enesis, the wide range of r e ported in vitr o v alues that can elicit EC r esponses suggest mor e gran ularity is needed to pinpoint specific ang iog enic behavior observed in vi v o in r esponse to specific hemodynamic forces.This is supported by other in vitro works that have shown strong connections between EC responses and TAWSS spatial gr adients [61][62][63][64] (ie , as opposed to single TAWSS v alues).Howev er, gi v en the difficulties of making such detailed measurements during ang iog enesis in vi v o, a critical first step in connecting in vi v o EC r esponses to spatial gradients is to dir ectl y identify 3D TAWSS spatial patterns that can exist during ang iog enesis, and hence a major focus of this work.The highr esolution sim ulations her e hav e ena b led quantification of these spatial variations at a highly localized level of detail.In terms of magnitude of variations, we have shown that the per-vessel TAWSS standard deviation can reach upwards of 55 dyne/cm 2 .Such v ariations ar e significant, gi v en the afor ementioned works, which have shown single TAWSS values themselves in this range to be physiologically relevant, let alone the spatial variations.
In addition to the magnitude of TAWSS spatial variations, we also quantified detailed maps of the TAWSSG directional components.The rich patterns further r ev eal and quantify the extent to which 3D spatial variations likely occur during ang iog enesis in vi v o.The existence of 3D spatial gradients over undulating surfaces of EC sheets was demonstrated in vitro by Barbee,7 who show ed gr adient hot spots in the flow direction can reach up to 7 d yne/cm 2 / μm.Ho w ECs respond to suc h gr adients has been r e ported by in vitro works, 61 , 62 which sho wed ho w EC migration patterns and directions are connected to flow orientation.In vi v o work by Franco et al. 65 has also shown EC migration directional behavior to be influenced by flow direction.Such directional response of ECs emphasizes the importance of elucidating the directional components of 3D TAWSSGs, which can occur in real ang iog enic vessel morphologies.Here we have shown that the TAWSSGs av era ged ov er each v essel can r each 15 dyne/cm 2 / μm (axial) and 13 dyne/cm 2 / μm (circumferential), with a significant heter ogeneity observ ed acr oss all v essels.Further considering the 3D TAWSSG spatial patterns w e observe , as w ell as maxim um v alues that can reach 300 dyne/cm 2 / μm at highly localized surfaces within vessels (see Supplementary Material), a new picture is painted of the directional forces (and their magnitudes) likely experienced by ECs in vivo during ang iog enesis.
Experiments showing connections between EC directional migration patterns and shear stress gradients have focused on flow direction gradients (ie, our axial direction).Yet, by modeling 3D vessels and RBC-influenced fluid dynamics our findings quantify significant circumferential gradients that can exist that are perpendicular to the flow direction.Earlier in vi v o works in non-ang iog enic micr ov asculatur es, such as those by Kim and Sarelius, 90 and Noren et al. 94 showed how shear stresses can be different on opposite sides of vessels, which suggests the existence of circumferential TAWSS variations.This was especially noted in the latter of these works within bifurcations.We have shown and quantified here the existence of such spatial patterns in 3D, and not only within bifurcations but also along capillaries due to local surface undulations arising fr om v essel cr osssectional variations and centerline curvature changes.Importantl y, we hav e shown that RBCs significantly enhance gradients in the ang iog enic networks yielding values roughly 3-5 times greater than without considering RBCs, and that the circumferential component is increased by 5-10 times compared to the axial component.These findings emphasize the role played by 3D geometric complexity and individual RBCs in driving the complex force profiles experienced by ECs during ang iog enesis in vi v o.
Micr ov ascular loop structur es ar e also shown to experience regions of highly localized TAWSS spatial variations.We observe loops to typically have at least one vessel with negligible TAWSS v alues, naturall y leading to sharp gradients at locations where such vessels connect to other vessels.Furthermore , w e have sho wn ho w some loop structur es can hav e v essels with m ultidirectional flow due to RBCs.Combined with patterns due to local 3D surface complexities, our findings generall y r ev eal loop structures to be regions with significant TAWSS spatial variations.Recently, Ghaffari et al. 38 reconstructed isolated micr ov ascular loops in 2D based on in vi v o ima ges fr om the capillar y plexus of avian embryos, and predicted TAWSS values ranging fr om negligib le to 1.4 dyne/cm 2 using a single-phase blood flow model.Nota b l y, they observ ed the locations of TAWSS minima, excluding regions where flow streams merge, to generally correlate with locations of new sprout growth.While TAWSS variations in loop structures predicted in the current work are nearly two orders of magnitude larger, and despite differences in the nature of the ang iog enic processes, these previous findings demonstrate the ability to connect TAWSS patterns to specific and highly localized ang iog enic events in vivo.Thus, the new 3D beha vior re vealed here encourages de veloping ne w mechanistic understanding of EC movement during ang iog enesis in vi v o based on high-resolution 3D surfaces and RBC-r esolv ed fluid dynamics.
Recent findings connecting sprout formation to regions of low TAWSS speak to the importance of identifying the existence of localized shear stress extremum in ang iog enic micr ov ascular networks along with influencing factors.We have shown here that 3D microvascular surface features and the presence of RBCs enhance both TAWSS hot spots and cold spots.Enhancement of cold spots by RBCs resulting in focused re gions of low ered TAWSS ( F igure 8 ) is a mechanism supporting our observation of RBC-enhanced TAWSS minima per vessel ( Figure 6 A-ii).TAWSS hot spots due to changing vessel curvatur e ( Figur e 8 B) r esult in distinct patterns alternating between vessel sides, with lower TAWSS values occurring on the apices of tortuous v essels.Inter estingl y the latter is in a gr eement with recent in vivo work that connected in vivo vessel tortuosity to incr eased spr outing on the apical sides. 71Hot spot enhancement is also observed here due to the wrapping motion of RBCs flowing around pillar-like obstructions ( Figure 8 A), a geometry that bears resemblance to that reported by Arpino et al. in an in vi v o study identifying such structures as inte gr al to intussusce pti v e ang iog enesis.The imag es of our ang iog enic micr ov ascular networks come from single points in time, and thus it is unclear whether the vessels in this region formed thr ough intussusce ption.Nev ertheless, our findings suggest asymmetrical shear stress patterns significantly enhanced by RBCs may influence behavior in ang iog enic g eometries of this nature.
Lastly, we have augmented the concept of WSS variations and heterogeneity by also analyzing variations in time.Examples of instantaneous WSS patterns corresponding to local RBC configurations and near-wall behavior suggest a causal role of RBCs on hot spots.This results in spatial patterns that move and ev olv e with the RBCs, and fluctuations in WSS statistics occur on timescales similar to that of individual RBC passage.To provide context for this, the overall physical time simulated is on the order of 1 s, during which time many RBCs (on the order of 1000) may pass by a particular wall location.This justifies our assumption of rigid vessel walls, considering the timescales of vessel adaptation are orders of magnitude longer.The same argument, however, suggests that the influence of the temporal variations r e ported her e on ang iog enic chang es may be questiona b le.Nevertheless, we include the time-dependent behavior to provide a thorough description of the model predictions and disseminate findings.
Overall, our findings reveal a new high-resolution 3D picture of the micr oenvir onment likel y underl ying ang iog enic processes in vi v o.This adv ancement pr ov okes new questions to understand roles of 3D shear stress patterning on EC behavior along vessels.Such provocation motivates future studies correlating the presence of 3D shear stress spatial and temporal patterns with patterns of EC phenomena, such as local hot spots and EC permeability or junctional disruptions, or spatio-temporal variations and signaling coordination.

Limitations and Additional Considerations
Appreciation of the presented findings requires discussions on the modeling limitations.The many big-picture-type challenges associated with getting the model "right" pr ov oke consideration of the model's purpose.The challenges and assumptions with this type of modeling include aspects such as selection of boundary conditions, initial placement of RBCs and injection at inlets to maintain hematocrit levels, resolving 3D geometrical features, or necessary compromises to manage long run times.A model is only as good as what it is used for, and for the current work, our goal is to identify potential shear stress patterns in remodeled network regions post angiogenesis taking into account features such as 3D fluid dynamics, RBC deformation and interactions, v essel connecti vity, and vessel specific curv atur es.The discussions below emphasize many of the c hallenges in ac hieving this, how w e hav e ov ercome them and resulting limitations, and provide context for interpretation of presented findings.
Boundary conditions for the simulations are prescribed flow rates at the inlet and outlet vessels.While we do not know what the actual boundary conditions should be, the flow values were chosen to yield shear rates spanning a physiolog ical rang e, 81 as discussed in Materials and Methods.Similarly, hematocrit values at inlets were prescribed following physiological values r e ported for the microcirculation. 81 , 82Given the known heterogeneity of hemodynamic quantities in the microcirculation, parameter selections that are also physiological could ostensib l y r esult in incr eased/decr eased shear str ess ma gnitudes and patterns beyond what is r e ported her e.Yet, v alues wer e chosen to span an "av era g e" physiolog ical rang e based on literatur e v alues and thus r esults pr ovide a r e pr esentation in this context of what may occur during ang iog enic remodeling.
The choice of flow versus pressure boundary conditions was specifically made to facilitate shear stress comparisons between RBC-r esolv ed sim ulations and plasma-onl y sim ulations.Fixed pr essur e boundar y conditions would r esult in similar per-v essel WSS values between the RBC and plasma-only simulations due to resulting similar per-vessel pressure drops based on the fluid mechanics, as observed in prior work. 8On the other hand, flow boundary conditions have enabled isolating contributions of RBCs to shear stress behavior.Fixed pressure boundary conditions could be used to war d understanding flo w distribution patterns and connections to hydrodynamic resistance, although we limit our focus here to the WSS alone.We also note that, even though we have flow boundary conditions, when we change them (either by the 0.5x or 2x factors) the per-vessel TAWSS throughout networks does not just linearly scale in proportion to the boundary flows.This is due to the presence of the RBCs. Figure 7 B shows how, if you have just plasma and you double or halve the prescribed flows at the boundaries, the shear stresses in all of the vessels scale in proportion.In contrast, if you have RBCs this is not necessarily the case, especially in smaller-diameter vessels.
The initial shapes of RBCs injected at inlets are deformed and highly asymmetric and derive from separ ate str aight tube simulations resulting in random placement over inlet cross-sections.After initial RBC injection, interactions with other RBCs as well as the complex vessel geometries quickly dominate RBC motion and dynamics.Different inlet positions of each injected RBC would likely affect the specific time-dependent WSS patterns in v essels, howev er, it is not expected that such differences would alter any of the main findings presented.
The physical dimensions of the modeled RBCs are consistent with human RBCs, 74 while the ang iog enic micr ov ascular networks consider ed ar e fr om the rat mesentery.The ultimate goal of this work is to understand and pr edict r e pr esentati v e shear stresses during ang iog enesis, and what we really care about is the human scenario.To this end, we modeled human RBCs, and the rat mesentery is justified because it displays characteristics of ang iog enic networks regardless of species.In terms of particle-scale effects on shear stresses, we note that human and rat RBC dimensions differ in their major radius by appr oximatel y 0.4 μm and that human RBCs have ∼30% greater surface area and 50% greater volume. 74 , 75This potentially would result in a greater impact of human RBCs on shear stress than rat RBCs.Further, based on these size differences, the number of rat RBCs r equir ed for a specific hematocrit level would be increased over human, but not by an order of magnitude.While different specific time-dependent WSS patterns would likely occur, the qualitati v e details of the r e pr esentati v e patterns shown here would not be fundamentally different.
The time-dependent WSS behavior reported here occurs entir el y because of the deforma b le RBCs, as the boundary conditions are fixed.This results in relatively high-frequency WSS oscillations ( ∼0.01 s) compared to larger-scale oscillations ( ∼1-100 s) due to various changes in vessel diameters attributed to local smooth muscle cell r eacti vity along arterioles, precapillary sphincters, 95 , 96 or ev en pericyte r egulation along capillaries. 97Additionally, work by Kiani et al. 98 both computationall y and experimentall y sho wed ho w RBCs themselv es can pr oduce oscillations in vascular networks with timescales in the range of ∼5-25 s.This warrants future consideration of RBCr esolv ed sim ulations at these longer timescales to investigate resulting oscillations and effects on WSS.The three different boundary flow strengths considered in the current work provide some measure of varying WSS responses to flow conditions, but the dynamics of tr ansitioning betw een flow states is not known.While the current w ork isolates tempor al fluctuations due to RBCs, interrelated effects of fluctuations due to other sources, both at vessel and sub-vessel levels, remain to be investigated.We also note that for each of the three flow conditions considered for eac h netw ork, the geometries remain fixed.For such scenarios, the predicted per-vessel values across these conditions ar e alter ed such that both the av era ge shear str ess v alues and the heterogeneity increase with flow strength ( Figure 7 A).In addition, the distribution of shear stress values among the v essels is alter ed in a non-linear fashion due to the presence of RBCs ( Figure 7 B).For other scenarios in which the simulated network diameters may change, this would likely alter the heterogeneity as well, although this remains to be investigated as well.
Due to the computational effort inv olv ed, sim ulations span a physical time on the order of 1 s.Effects of oscillations at the longer timescales noted a bov e ar e difficult to estimate.On one hand, due to the flow regime it might be justified to assume the behavior can be captured by a superposition of separate sim ulation r esults at differ ent, steady boundar y conditions.Yet on the other hand, due to the presence of RBCs, and considering changes in micr ov ascular tone/structur e ov er time, ther e are non-linear effects present that could influence the behavior during micr ov ascular transitions.An important point to note though is that, in a broad sense, what this work is highlighting is localized 3D heterogeneity.If one were to consider a 1 s window 10 min later, the specific results might be different, but this heterogeneity will still exist.The real extent of the spatiotemporal heterogeneity though, is unknown.So, while we do not know exactly if or how such oscillations will affect behavior, future inv estigations ar e w arranted.
Vessel tortuosity is a hallmark ang iog enic micr ov ascular c har acteristic, and here we have shown how shear stress hot spots form in tortuous vessels near regions of high surface curv atur e ( Figur e 8 B), similar to other r ecent work. 8Specific c har acteristics of tortuous vessels are expected to vary based on location in the body and function.The findings here provide evidence that the ang iog enic tortuous vessel c har acteristic results in distinct 3D shear stress patterns and warrants future investigation to identify tissue and organ-specific beha vior.F or example, in skeletal muscle with increased vessel tortuosity in conjunction with rapidly changing RBC fluxes, there may be unique enhancements to shear patterns to achieve a functional end, although future investigations are needed.
For the current study, the selection of ang iog enesis stimulation w as arbitrar y; gi v en the study's focus on modeling RBC flow through ang iog enic netw orks, r ationale for the use of compound 48/80 stimulated tissues is supported because the networks display ang iog enic remodeling c har acteristics.While the specific mechanism remains unknown, the stimulation has been associated with mast cell release of histamine and other highly acti v e substances. 84 , 85It is unknown whether WSS is a triggering mec hanism.How ever, based on these previous studies, we speculate that ang iog enesis is initially triggered by changes in the local chemical environment.An intended main takea wa y of the current study is the identification of potential RBC flo w d ynamics, including local WSS spatiotemporal patterns and the use of the acti v el y r emodeling network r egions ena b les this.The r e pr esentati v e r egions used for the curr ent study wer e selected because they displayed high vessel density and patterns r eflecti v e of acti v e network r emodeling.Regions of interest did not include capillary sprouts.Hence, the vessel patterns in the ang iog enic reg ions can be assumed to be post the initial capillar y spr outing sta ge.Futur e studies will be needed to evaluate RBC flow effects due to sprout presence, sprout dynamics, and sprout connection with other vessels.The lack of sprouts in the regions for the current studies emphasizes that our r esults pr ovide insights into the effect of remodeling on local RBC flow dynamics versus the implication of RBC flows on sprouting.Further, our results motivate similar future studies to determine whether similar RBC flow dynamics are associated with ang iog enic patterns due to other alternati v e stim uli and acr oss differ ent tissue beds.
This work focuses on the shear stresses experienced due to hemodynamics and does not include other factors known to contribute to ang iog enesis, such as tissue growth factors, interstitial transport, or transendothelial transport.Following the discussion a bov e r elated to time-de pendent behavior, the timescales r esolv ed her e ar e of a different order than those associated with these processes, and thus impacts on our findings are expected to be small.In terms of 3D vessel surface features, the model also does not account for the endothelial surface layer (ESL) or individual EC surface structures.3D networks were constructed assuming locally circular lumen shapes and do not r esolv e asymmetrical cr oss-section shapes or surface irregularities.While such considerations will likely influence the exact 3D shear stress patterns, here we have focused on revealing types of patterns and behavior that can occur when 3D ang iog enic vessel structures and RBCs are resolved.
The findings suggest that including individual EC geometries or other ang iog enic surface irregularities may lead to unique 3D shear stress patterns.For example, vessels with changing crosssectional shapes and out-of-plane asymmetries may result in an increased frequency of local hot and cold spots, especially wher e mor e sev er e local surface v ariations occur.This is supported by the patterns we have shown in Figure 4 B near the red and blue arrows, which are similar in spirit to this.Additionally, inclusion of individual EC shapes may result in sharp WSS gradients at locations where ECs connect to one another, due to locally high surface curvatures near EC junctions.In the current work, high curv atur e of local surfaces such as Figure 8 B (feeding v essel), r esulted distinct hot spots and large local surface gradients.Going forward, the specific details and idiosyncrasies of such patterns are difficult to predict without explicitly modeling them.We have shown here that revealing such details requires high-r esolution RBC-r esolv ed sim ulations, yet in the future, the output from such simulations could in principle be used to build impr ov ed r educed-order models and circumvent the need for large-scale simulations.Altogether, the identification of the r e pr esentati v e patterns presented here is an important first step needed to inspire new resear c h directions for vascular biologists and modelers alike to tackle ideas such as ESL/EC shape influences.
and reasoning behind the use of human RBCs is provided in the Limitations and Additional Considerations section.The governing fluid flow equations are the unsteady Stokes equations for a constant density ( ρ) varia b le viscosity ( μ) fluid ρ

Figure 2 .
Figure 2. 3D RBC-r esolv ed sim ulations in ang iog enic micr ov ascular networks.(A) Ima ge data for network 1, scale bar is 40 μm; (B) snapshot of RBC-r esolv ed sim ulation; (C) surface contours of the time-av era ged WSS, (D-F) same information but for network 2, (G) predicted time-averaged WSS (TAWSS) per vessel versus diameter, (H) heatmap of TAWSS per vessel as a function of vessel length and diameter.

Figure 3 .
Figure 3. Vessel to vessel T AWSS heterogeneity .Representative images of 3D TAWSS surface contours predicted by simulations, illustrating the large degree of TAWSS v ariation fr om one v essel to another.(A) High TAWSS capillaries (r ed arr o ws) connecting to lo w TAWSS v en ule; (B) gr oups of low TAWSS v essels (b lue arr ows) connected to groups of high TAWSS vessels (red arrows); (C, D) high and low TAWSS vessels interconnected.

Figure 4 .
Figure 4. Ang iog enic TAWSS spatial variations and gradients in 3D.(A) TAWSS contours for a r e pr esentati v e ang iog enic netw ork re gion, scale bar is 10 μm, and arr ows gi v e the flow dir ection; (B) zoomed in view of a subr egion fr om (A), detailing the complex TAWSS spatial patterns and significant v ariations, arr ows point to r e pr esentati v e hot and cold spots, scale bar is 5 μm; (C) standard deviation ( σ ) of TAWSS versus vessel diameter, where each data point gives σ for all surface points in a vessel; (D) spatial gradients of the TAWSS components ( ∇ ax , ∇ circ ) v ersus v essel diameter.Each data point gi v es the r especti v e ∇ component av era ged ov er all surface points in a vessel.(E, F) Contours of the TAWSS spatial gradients (E, axial direction; F, circumferential direction) corresponding to the TAWSS patterns in (B).

Figure 5 .
Figure 5. Spatial patterns around microvascular loops and average TAWSS distribution behavior in vessels.(A) Representative examples of microvascular loop structures and resulting TAWSS spatial patterns.Black arrows give the flow direction, yellow arrows and the X symbol denote a low flow m ulti-dir ectional v essel and a no-flow v essel, r especti v el y.Gr een arr ows point to transition regions between connected vessels with significant TAWSS variations over short distances.Scale bar is 6 μm.(B) Histogram quantifying tendency of vessels to experience low or high TAWSS, on a per-area basis.Data are averaged over all vessels.

Figure 6 .
Figure 6.Alteration of TAWSS behavior in ang iog enic networks due to RBCs.(A) Ratio of TAWSS per vessel from simulations with and without RBCs ( TAWS S RBC / TAWS S plasma ) v ersus v essel diameter.(i) Ratio based on mean, maximum, and standard deviation of TAWSS per vessel, (ii) ratio based on minimum TAWSS per vessel, with values below 1.0 indicating vessels where the minimum TAWSS decreases due to RBCs; (B) TAWSSG per vessel with and without RBCs, versus vessel diameter.(i) Ratio α = TAWSS G RBC / TAWSS G plasma for both circumferential and axial components, (ii) ratio of α circ /α axial versus diameter.

Figure 7 .
Figure 7. Influence of boundary conditions on TAWSS and influence RBCs.(A) Distribution of TAWSS values across networks for the baseline flow conditions at network boundaries (BASE, or γ = 1.0), as well as values for flows modified from the base by factors of 0.5 and 2.0 ( γ = 0.5,2.0).(B) Ratio of TAWSS per vessel for each γ to the baseline conditions versus diameter.Left-and right-most plots give data from simulations with and without RBCs, respectively.

Figure 8 .
Figure 8. TAWSS hot and cold spots and enhancement by RBCs.(A) TAWSS field for a r e pr esentati v e r egion, with specific hot and cold spots discussed in the main text identified by the red and blue arrows, respectively.Also shown is a snapshot of RBCs, depicting representative behavior that gi v es rise to the noted hot and cold spots.(B) Hot and cold spots within a r e pr esentati v e bifurcation fr om sim ulations with and without RBCs The b lue arr ow points to an enhanced cold spot, which has a lower av era ge v alue and is mor e focused with RBCs than without.Red arrows identify hot spots due to surface curv atur e changes, with hot spot intensity enhanced by RBCs.

Figure 9 .
Figure 9. Time-dependent WSS behavior.(A) Snapshot of RBCs at an instant in time.Inset frames track the motion of the circled RBC as it moves through the vessel, with corresponding instantaneous WSS patterns.(B-D) Time signals for WSS statistics from two representative vessels circled in (A).(B) Mean vessel WSS, (C) maximum vessel WSS, and (D) standard deviation and minimum vessel WSS.r e pr esentati v e examples that depict and quantify aspects of the time dependent WSS predicted by the 3D simulations.To illustrate the relationship between instantaneous WSS spatial patterns and the movement of individual RBCs, the inset to Figure9A depicts WSS contours at three sequential time points along with the corresponding RBC positions.At each time point, the WSS contours exhibit spatial variations, and these can be observed to change following the RBC motion.This observation augments the concept of WSS heterogeneity by revealing how the sub-EC scale spatial variations are constantly in flux.To quantify the degree of WSS variations in both space and time , F igur e 9 B-D pr ovide example time signals for different WSS spatial statistics from two separate vessels identified in Figure9A.To re veal beha vior resulting from different hemodynamics expected in ang iog enic networks, we focus on (1) a small capillar y wher e RBCs flo w in single file and (2) a lar ger diameter v essel wher e RBCs flow in m ultifile fashion.The figur es focus on behavior over a time interval of 0.07 s in order to highlight the timescale of the fluctuations in quantities.First, shown in Figure9B is the time signal for the mean WSS (ie, WSS spatiall y av era ged ov er all points in the vessel at each time instant).As can be seen, the small capillary fluctuates r oughl y between 45 and 55 dyne/cm 2 , with an oscillation timescale of appr oximatel y 0.021 s.We calculate this timescale as the inverse of the dominant frequency resulting

Figure 10 .
Figure 10.Time-dependent WSS fluctuations within vessels.RMS values of the time dependent fluctuations of WSS statistics versus vessel diameter.The left-most plot gi v es RMS v alues for the mean and standard deviation WSS time signals for each v essel, while the right-most plot gi v es RMS v alues for the maxim um and minim um WSS signals.Best-fit lines corresponding to each data set are also provided based on a power-law r egr ession.