Association between erythrocyte dynamics and vessel remodelling in developmental vascular networks

Sprouting angiogenesis is an essential vascularization mechanism consisting of sprouting and remodelling. The remodelling phase is driven by rearrangements of endothelial cells (ECs) within the post-sprouting vascular plexus. Prior work has uncovered how ECs polarize and migrate in response to flow-induced wall shear stress (WSS). However, the question of how the presence of erythrocytes (widely known as red blood cells (RBCs)) and their impact on haemodynamics affect vascular remodelling remains unanswered. Here, we devise a computational framework to model cellular blood flow in developmental mouse retina. We demonstrate a previously unreported highly heterogeneous distribution of RBCs in primitive vasculature. Furthermore, we report a strong association between vessel regression and RBC hypoperfusion, and identify plasma skimming as the driving mechanism. Live imaging in a developmental zebrafish model confirms this association. Taken together, our results indicate that RBC dynamics are fundamental to establishing the regional WSS differences driving vascular remodelling via their ability to modulate effective viscosity.


Introduction
Sprouting angiogenesis is an essential vascularization mechanism and consists of two well-differentiated phases: sprouting and remodelling [1,2]. During the sprouting phase, a primitive network of vessels is laid out in response to proangiogenic factors via a well-established programme of cellular and molecular events [3]. The remodelling phase is responsible for overhauling this primitive network into a hierarchical structure that can efficiently implement the transport capabilities of the cardiovascular system. During the remodelling phase, extensive vessel pruning is achieved primarily via dynamic rearrangement of endothelial cells (ECs) [4].
In previous works, we demonstrated that blood shear stress coordinates EC migratory behaviour to achieve vessel regression during the remodelling phase [4,14,18]. In particular, differences in wall shear stress (WSS) between neighbouring vessel segments lead to polarization and migration of ECs away from vessel segments experiencing low WSS. In these studies, WSS was calculated using a mathematical flow model that assumes blood as a single-phase non-Newtonian fluid (i.e. homogeneous fluid of variable viscosity) without considering the presence of individual red blood cells (RBCs). However, recent computational studies in microscale vessels have demonstrated that RBCs leave transient WSS luminal footprints and therefore could non-trivially modify the local WSS differences driving vascular remodelling [23][24][25]. This effect of RBCs is closely related to their crucial role in the formation of cell-free layer and the regulation of effective viscosity as well as flow field in microvascular blood flow [26][27][28][29][30][31][32].
In the current study, we propose that the cellular nature of blood (i.e. primarily a suspension of RBCs) plays a key role in establishing the local WSS differences driving vascular remodelling. We approach the problem computationally based on the mouse retina model of angiogenesis and provide experimental validation in a developmental zebrafish model. We present simulations of cellular blood flow in vascular networks undergoing remodelling and characterize the bulk flow and RBC dynamics in them. Remarkably, we uncover a previously unreported high-level heterogeneity in RBC perfusion throughout the developing network and a strong association between RBC hypoperfusion and vessel regression, which we further confirm experimentally in our zebrafish model via live imaging. Furthermore, our experiments confirm previous findings with additional insights that the presence of RBCs is required for effective vascular remodelling at a whole plexus level [7]. Finally, we demonstrate that RBC hypoperfusion is primarily driven by the plasma skimming effect, i.e. the uneven distribution of RBC volume fraction at microvascular bifurcations [33,34], but uncover important deviations from existing theory caused by asymmetry of the cross-sectional distribution of haematocrit in some feeding vessels.
In line with our previous findings on WSS-modulated EC migration, we argue here that RBC hypoperfusion constitutes a mechanism for the enhancement of local WSS differences driving vascular remodelling. This is attributed to the direct relationship between RBC volume fraction, effective viscosity and WSS. Additionally, we speculate that vascular remodelling driven by the principle of removing RBC-poor vessels from the primitive vasculature will lead to a network layout that avoids portions of the tissue being vascularized but with poorly oxygenated blood. This RBC-driven process, highly dynamical and emergent in nature, can importantly contribute to the optimal patterning of vascular networks during development. Conversely, it provides a vascular remodelling mechanism to support recent clinical findings in diabetes mellitus and hypertension [35][36][37][38] reporting that the onset and progression of microangiopathy may arise from altered RBC mechanical properties (e.g. impaired deformability) that disrupts the microcirculatory blood flow: a hypothesis first proposed by Simpson in the 1980s [39].

Simulation of cellular blood flow in microvascular network: validation versus in vivo data
The vascular plexus of a wild-type mouse retina at postnatal day 5 (P5) was imaged and binarized following staining for collagen IV matrix sleeves (Col.IV) and ICAM2 luminal reporter (see protocol in §5.1.1). The ICAM2 mask delineates the perfused vessels in the network while the pixel-wise difference between the Col.IV and ICAM2 masks highlights vessel segments undergoing remodelling (figure 1a) as demonstrated in Franco et al. [4]. Therefore, the Col.IV mask constitutes a good approximation of the network morphology prior to these remodelling events. Indeed, the network recapitulates the lognormal distribution of vessel diameters previously reported by Bernabeu et al. [43], with a maximum value of 45 µm and a mean of 11.85 µm (figure 1b). Based on our previously proposed approach [44], we construct a threedimensional (3D) flow model from the Col.IV binary image and run whole-plexus blood flow simulations under the assumption of generalized non-Newtonian blood rheology (see the methods in §5.2.2). This is followed by simulations of RBC flow in selected regions of interest (ROIs; see the selection criteria in §5.2.1) from the capillary bed (e.g. ROI-1, ROI-2, ROI-3 and ROI-4 in insets of figure 1a) with inflow/outflow boundary conditions specified from the whole plexus model (see the methods in §5.2.3). For full details of the computational framework including model configuration and simulation parameters, refer to figures S1-S2 and tables S1-S3 in §S1 of the electronic supplementary material.
Our simulation recapitulates flow rates in the main artery and veins (running in the axial direction of the retina) of 0.36 µl min −1 and 0.17-0.19 µl min −1 , respectively. These flow rates are in good agreement with in vivo measurements recently performed in adult mouse retina of 0.39-0.59 µl min −1 and 0.24 µL min −1 , respectively [41]. The RBC velocities and fluxes calculated in the simulated ROIs are in the range 0-12.5 mm s −1 (figure 1c) and 0-1500 cells s −1 (figure 1e), respectively. These results are also validated against reported in vivo single-cell velocimetry data available for capillaries within diameter 3-7 µm [41,42]. Both the experimental and simulation RBC velocities follow a skewed distribution and show good agreement with median values of 2.14 mm s −1 and 2.35 mm s −1 , respectively (p > 0.05, figure 1d). The RBC fluxes are significantly higher in experiments than in simulations with median values of 120 cells s −1 and 19.3 cells s −1 , respectively (p < 0.01, figure 1f). The main reason for the lower median observed in simulations is the presence of a number of vessel segments with zero or negligible RBC fluxes in the primitive vasculature of our developing retina (figure 1f), which were absent in the reported experimental data for adult mouse retinal capillaries by [41,42]. This major difference will become the focus of study in the following §2.2.

Strong association between RBC hypoperfusion and vessel regression in developing vasculature of mouse retina
Having provided evidence of the broad range of RBC fluxes in the selected ROIs, we now investigate a potential association between RBC perfusion (referring to the state of vessels being either sufficiently or scarcely perfused by RBCs) and vessel regression within the network. We classify vessels in each ROI into three groups: lumenized, stenosis and regression (figure 2a  (d,f ) Comparison of simulated RBC velocities and fluxes against recent in vivo measurements by [41,42] from adult mouse retinal capillaries with the same diameter range (two-sided Mann-Whitney U test).
negative signal in ICAM2 for the vessel segment, indicating partial or complete lumen collapse.
In the four ROIs extracted from the whole plexus (including 63 vessels for ROI-1, 23 vessels for ROI-2, 28 vessels for ROI-3 and 30 vessels for ROI-4), we identified 27 regressed vessels, accounting for 28% out of 97 regressed vessels in total throughout the remodelling region of the plexus (see figure 1a). If the regressed vessels and vessels undergoing stenosis within each ROI are counted as remodelling events, the percentages of vessel remodelling within each ROI are 30.4%, 36.7%, 35.7% and 30.2%, respectively. The consistent level of remodelling events within the selected ROIs reflects that the overall remodelling of the primitive vasculature is relatively uniform from region to region.
Further examination of our cellular flow simulations reveals that most poorly RBC-perfused channels (hereafter defined as 'RBC hypoperfusion') in the simulation coincide with vessels in the regression group ( figure 2a,b). To quantify this, we record the trajectories of all RBCs within each ROI throughout the simulation and study their density across the ROIs (figure 2c and figure S5 in the electronic supplementary material). In general, a vessel segment with higher density represents good RBC perfusion, whereas those with low density indicate RBC hypoperfusion. Subsequently, the timeaverage RBC flux within each vessel segment is calculated and assigned to the three groups under study (figure 2d-g). Our analysis demonstrates that vessel segments in the lumenized group of each ROI have significantly higher RBC fluxes than those in the regression group (p < 0.05 for all ROIs, figure 2d-g). Meanwhile, the difference in RBC flux between the stenosis group and the regression group is not significant (p > 0.05) for 3 out of 4 ROIs (except for ROI-1 due to a single outlier, figure 2d-g). Furthermore, there is no significant difference in the distribution of RBC fluxes between the lumenized groups of these ROIs, despite different maximum values recorded depending on the relative location

In vivo validation of the effect of RBC perfusion on vascular remodelling in zebrafish caudal vein plexus
To provide experimental confirmation of the association between RBC hypoperfusion and vessel regression predicted by our computational model, we turned to a zebrafish model of vascular development, where simultaneous live imaging of vessel remodelling and RBC dynamics is possible. We chose the caudal vein plexus (CVP) for observation 48-72 h post fertilization (hpf), a period during which gradual remodelling of the plexus down to a single, well-defined vascular tube begins [45]. Comparison was made between control (ctl) morpholino oligomer (MO) fish with normal RBC perfusion (figure 3a) and gata1 MO fish not carrying RBCs in the bloodstream (figure 3b, see §5.1.2 for experimental details). Our time-lapse imaging of the CVP in ctl MO fish captures heterogeneous RBC perfusion (figure 3c) leading to multiple findings of intermittent and complete RBC depletion in vessel segments (figure 3d), followed by vessel stenosis (figure 3e) and eventual regression (figure 3f). These in vivo findings therefore confirm our computational predictions. Next, we investigated how RBC deletion in gata1 MO fish impacts CVP remodelling at a network level. We measured CVP widths at standardized locations along the anterior-posterior fish axis given by the positions of eight consecutive intersegmental vessels (ISVs), for both the ctl MO and gata1 MO groups (seven fishes in each group) at 50 hpf and 72 hpf (see figure 3g,h for an example). During this period of time, substantial remodelling leading to narrowing of the plexus is observed in the wild type (agreeing with [45]). Furthermore, significantly larger reduction of CVP width in the ctl MO   royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210113 group is found in comparison to the gata1 MO group, accounting for 9% and 1% of relative width change (mean reductions 7.64 µm and 2.01 µm), respectively (p < 0.01, figure 3i). To exclude that differences in fish growth may confound these results, the longitudinal extension of the somites is also examined, for which no significant difference is found between the ctl MO and gata1 MO groups (p > 0.05, figure 3j). Taken together, this several-fold difference in CVP width change implies that the presence of RBCs is necessary for normal CVP remodelling 50-72 hpf, and therefore the heterogeneity in RBC perfusion described in figure 3c-f does play a role in orchestrating network-level remodelling.

RBC depletion in the network not predictive by vessel diameter
The absence of RBCs in some vessel segments in both our simulations and experiments poses a crucial question about cellular flow in developmental vascular network: what is the governing mechanism that determines which vessels to be perfused with cells and which to be devoid of? It is tempting to speculate that the pattern of vessel RBC hypoperfusion is merely a sizeexclusion effect; namely, certain vessel segments are simply too narrow to allow cells to pass through. However, the vessel diameters encountered in the present simulations (about 2-20 µm, see figure 1c,e) are unlikely to be the dominant factor affecting RBC transit, since the high deformability of RBCs under physiological conditions enables them to pass through exceedingly small passages as narrow as 1-2 µm [46][47][48].
Indeed, the material model we adopt for the RBC membrane in simulations (see benchmark tests in figures S3-S4 and §S1.3 of the electronic supplementary material) allows the cell to adopt highly elongated shapes to pass through nearly all narrow capillaries that they enter (with two occasions of cell occlusion due to negligible flow in the vessel), whereas some larger vessels are devoid of RBCs (see the electronic supplementary material, movies S1-S4).
To further quantify the pattern of RBC perfusion in the network, we define RBC 'depletion' and 'enrichment' by examining the sign of the term DQ Ã ¼ Q Ã rbc À Q Ã blood (see definitions in §5.3.1) for child branches at diverging bifurcations, which evaluates the fractional RBC flux a child branch receives from its parent Q Ã rbc relative to its counterpart of fractional blood flow Q Ã blood . If ΔQ* is negative for a branch, it is regarded as 'RBC-depleted' as disproportionately fewer RBCs are allocated to it than the proportion of blood flow it receives; otherwise the branch is defined as 'RBC-enriched'. By inspecting ΔQ* against vessel diameter D vessel for all valid child branches extracted from ROI-2, ROI-3 and ROI-4 (44 vessels in total, electronic supplementary material, figure S5), we find that RBC depletion occurs throughout the whole range of vessel sizes investigated (D vessel ∈ [2, 16] µm, figure 4a). Notably, for one child branch with D vessel ≈ 9 µm (larger than the physiological RBC diameter), ΔQ* equals −0.2, indicating a 20% reduction of RBC transit relative to the blood flow in it. Based on these findings, we conclude that RBC depletion is not a size-exclusion effect. Meanwhile, it is found that RBC enrichment happens only in medium/large vessels (D vessel > 5 µm). Within the intermediate diameter range D vessel ∈ [5,12] µm, the vessels have nearly equal chances of being enriched or depleted.

Plasma skimming as a mechanism for RBC depletion in developing vascular network
Having demonstrated that size-exclusion effect is not sufficient to explain the RBC depletion observed in §2.2, we now turn our attention to potential haemodynamic mechanisms. First, we inspect the RBC fluxes in the vessels against their blood flow rates, and find the two variables nearly monotonically correlated (see figure S6b in the electronic supplementary material). Furthermore, we observe that characterizing DQ Ã ¼ Q Ã rbc À Q Ã blood against the haemodynamic indicator Q Ã blood satisfactorily separates the RBC-depletion zone from the RBC-enrichment zone (figure 4b). This finding is in line with the plasma skimming effect described by empirical models such as the widely employed phase-separation model (PSM, see introduction in §5.3.2), first proposed by Pries and co-workers [49,50].
To qualitatively assess the agreement between our data and the PSM, we plot the fractional RBC fluxes Q Ã rbc in individual child branches (44 studied here) of any divergent bifurcation against their fractional blood flow Q Ã blood (figure 5a). The distribution is indeed reminiscent of the sigmoidal relationship predicted by the PSM, where Q Ã rbc . Q Ã blood for Q Ã blood . 0:5 up to a threshold where Q Ã rbc ¼ 1 (and the opposite for Q Ã blood , 0:5 down to Q Ã rbc ¼ 0). Furthermore, albeit with occasional exceptions, the child branch relatively smaller in  size (red circles) within a bifurcation tends to receive lower blood flow and consequently fewer RBCs, whereas the larger child branch (blue squares) is more likely to have higher blood flow and attract more RBCs (inset of figure 5a). Note that 'smaller' or 'larger' here is a relative notation between the two child branches within a bifurcation, instead of a measure of the absolute vessel size. Furthermore, we compare our simulation data at each bifurcation with empirical predictions given by the PSM (equation (5.1) and electronic supplementary material, equations (S1)-(S3)). We observe less than 5% errors for 18 out of 22 bifurcations (see figures S7-S8 for individual comparison and table S4 for complete error evaluation in §S2 of the electronic supplementary material). Quantitatively similar results are obtained if the potential effect of absolute cell volume difference between the RBC model and realistic mouse RBC is also accounted for (see formulation in equations (S4)-(S6), and complete error evaluation against simulation data in table S5 of the electronic supplementary material). Figure 5b-e demonstrates four cases where good agreement is shown. In the first bifurcation, both child branches have considerable proportions of blood flow (roughly 30% and 70%, respectively) and are well perfused by RBCs, with the fractional RBC fluxes matching the PSM predictions (figure 5b). In the second bifurcation, most RBCs (Q Ã rbc % 90%) enter the relatively larger child branch as it receives more than 80% of the blood flow from the parent branch (figure 5c). In the third and fourth bifurcations (figure 5d,e), the smaller child branch is nearly devoid of cells as the relatively larger branch attracts almost all RBCs from the feeding vessel owing to its predominantly higher proportion of blood flow (Q Ã blood . 85%). Larger than 5% deviations from the PSM are observed in the remaining 4 of the 22 bifurcations studied. Shown in figure 6a,b are two cases where the simulated RBC fluxes in the bifurcation deviate from the empirically predicted values. Interestingly, the RBC flux contrast between the two child branches is slightly underestimated by the PSM in the first case (figure 6a) whereas it is substantially overestimated in the second case (figure 6b). To investigate these disagreements, we examine the flow streamlines in the midplane of the bifurcation and calculate the relative distance χ of the stagnation streamline separating the flow entering the left branch from that entering the right branch ( figure  6c,d ). Clearly, the higher-flow branch (child branch with larger Q Ã blood ) receives proportionally more streamlines from the parent vessel and should therefore receive more RBCs accordingly provided that the cells are axisymmetrically distributed in the parent branch (as assumed by the empirical PSM).
Next, we investigate whether the assumption of haematocrit axisymmetry by the PSM holds in our simulations. To this end, we calculate the cross-sectional distribution of RBCs in the parent branch. We find clearly asymmetric distributions, especially for case two (figure 6e,f ). In both cases, the haematocrit profile is skewed towards the right-hand side of the bifurcation, with a cumulative haematocrit to the left-hand side of the vessel centreline either slightly or substantially smaller than 0.5 (see insets of figure 6e,f ). Such haematocrit asymmetry makes the downstream child branch on the right-hand side inherently advantageous for RBC intake (hereafter referred to as haematocrit-favoured branch). In case one, the higher-flow branch coincides with the haematocrit-favoured branch, resulting in enhanced RBC flux difference between the two child branches; whereas in case two, the higher-flow branch differs from the haematocrit-favoured branch, thus attenuating the RBC flux difference. A full description of the mechanism leading to this skewness, which originates from the interplay between the complex geometry and the emerging RBC behaviour in the microvasculature, is out of the scope of the present study and will be explored in future work.

Discussion
Earlier studies have extensively explored the effect of blood flow on cardiovascular development and it is well accepted that haemodynamic cues are essential for vascular development [7,11,12]. Previous work by our groups further established the pivotal role of regional WSS differences, which were found to locally modulate the polarized migration of ECs into highshear vessel segments and cause the regression of adjacent segments experiencing low shear [4,14,18]. However, because the quantification of WSS in these studies relied on mathematical flow models assuming simplified blood rheology, the question of how the presence of RBCs and their profound impact on microvascular haemodynamics affect vascular patterning has not been addressed. Recent simulations of cellular blood flow in single microvessels have shown that RBCs can non-trivially modify both the mean and oscillatory components of local WSS [23,25], thereby suggesting a salient impact of RBCs on the mechanotransduction of fluid forces during angiogenesis. Based on these findings, we hypothesize that RBCs in the developing mouse retina play an active role in the course of vascular patterning towards a functional network via altering WSS spatially and temporally.
In the current study, we simulated RBC dynamics in the vascular plexus of a wild-type mouse retina at P5 ( postnatal day 5). The computed RBC velocities are in good agreement with in vivo measurements made in adult mouse retinas [41,42]. However, differences arise in terms of RBC fluxes (more heterogeneous distribution with lower median value, figure 1f ). We attribute this discrepancy to structural differences between the developmental stage simulated and the adult stage considered for validation. We hypothesize that as the primitive network remodels, the simulated RBC fluxes will become closer to the in vivo measurements from adult mice. Future advances in live imaging of the mouse retina model will contribute to elucidating this process.
There are other sources of uncertainty in the current simulations remaining to be quantified. First, the discharge haematocrit at the inlets of all ROIs was fixed at 20% for consistency. This value is in line with the average microvascular haematocrit level reported in the literature, i.e. 23 + 14% by [52], but may not necessarily reflect the haematocrit levels of individual mice studied in [41,42] and used for validation in our study. Second, the flow in realistic microvessels with smaller lumen size than the RBC diameter is typically  Figure 6. Occasional deviation of simulation data from the empirical model [51] due to the asymmetry of haematocrit profile in the parent branch. royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210113 impeded due to the close contact of RBCs with the endothelial surface layer (ESL, consisting of glycocalyx and/or cilia), whereas in our present in silico model such microscopic structures on the vessel wall are not considered. Prior studies in glass tubes where the ESL was absent demonstrated a twofold decrease in flow resistance when compared with in vivo conditions [53,54]. Our simulations of cellular blood flow reveal a previously unreported high-level heterogeneity in RBC perfusion throughout the developing vascular network, where a number of plasma vessels exist with rare RBC transit over time, and a strong association between RBC hypoperfusion and vessel regression. From a mechanistic point of view, the effective viscosity in neighbouring RBC-depleted and RBCenriched branches can differ substantially. Such disparity in viscosity will enhance the regional WSS difference between neighbouring branches, which can in turn promote the pruning of the vessel segment depleted of RBCs according to our previously reported mechanisms [4].
We provide further experimental confirmation of the above findings in a developmental zebrafish model, which is amenable to live imaging. In agreement with Lucitti et al. [7], we show that the presence of RBCs is necessary for effective remodelling at a whole plexus level. Furthermore, we extend the conceptual model by demonstrating that intermittent and complete RBC depletion (also observed recently in a model of inflammatory corneal angiogenesis [55]) selects vessels for pruning, which are likely to experience low shear and become unfavoured during vascular remodelling. Apart from their impact on the WSS differences, Xiong et al. [56] recently reported that the signalling role of RBCs also influences vessel remodelling as erythrocyte-specific sphingosine-1-phosphate was found crucial for vascular stabilization and maturation during embryonic development. Understanding the interplay between haemodynamic forces and secreted factors derived from RBCs, however, is beyond the scope of the present study.
To fully understand the perfusion of RBCs within the developing mouse retina, we asked the question: what is the mechanism behind RBC depletion/enrichment in the primitive network? Quantification of RBC perfusion against vessel diameter and blood flow in individual bifurcations rules out vessel-size exclusion as the primary factor. Instead, the empirical model by Pries et al. [49,51] based on the plasma skimming effect explains the uneven partitioning of RBC fluxes satisfactorily in 18 out of 22 cases, therefore implying that the distribution of RBCs within the developing network is flow-mediated rather than geometry-dominant.
For the four cases where predictions of the empirical model do not satisfactorily match the simulation data, we find considerable haematocrit asymmetry in the cross-sections of the feeding branches, which is against the central assumption of axisymmetric haematocrit profile by the model [49,51]. This makes accurate prediction of the RBC perfusion in a given vascular network challenging without certain knowledge of cross-sectional cell distributions. Our observation of haematocrit-favoured and haematocrit-unfavoured vessel branches at microvascular bifurcations is in line with recent in vitro findings of inversion of the classic haematocrit partitioning, which otherwise always favours the higher-flow branch [57][58][59][60]. Similar reverse partitioning was reported by simulations of cellular blood flow in microvascular networks with vessel diameters designed following Horton's Law [30]. Our group recently identified reduced interbifurcation distance and complex branching topology as sources of haematocrit bias in the context of tumour blood flow, where a role was proposed for the resulting abnormal partitioning in establishing tumour tissue hypoxia [61].
The asymmetry not only applies to the haematocrit profile but also the velocity profile. For a majority of microcirculatory models, the Poiseuille Law has been employed to simplify the haemodynamics and reduce computational cost, manifested by parabolic velocity profiles featuring a peak velocity at the centreline of microvessels. Our results counter such a simplification, as the velocity profile in a capillary vessel can significantly deviate from a parabola and become skewed over time in the presence of travelling RBCs (figure S9 in §S3 of the electronic supplementary material). Similar deviation of velocity profiles has also been confirmed by imaging data from living mouse retina, where a 39% error in flow estimation was reported if assuming a parabolic profile [41].
Finally, the implications of the uncovered association between RBC hypoperfusion and vessel regression are not limited to developmental vascular remodelling. Diseases such as diabetes mellitus and hypertension lead to changes in RBC deformability [62][63][64] and these have been associated with vascular complications such as diabetic retinopathy or nephropathy in several clinical studies [35][36][37][38]. However, the magnitude of the changes in RBC deformability reported in the literature (ranging between 10% and 50% [38,65,66] to several folds [67,68]) would not support a model where such biomechanical changes alone are sufficient to cause vessel occlusion since previous studies showed that order of magnitude changes in shear modulus are required to impede RBC transit through narrow passages [46]. By contrast, our findings support a new concept where changes in the mechanical properties of the RBC membrane leading to abnormal haematocrit partitioning at bifurcations (via altered radial distributions of RBCs [69][70][71]) would reintroduce, in adult networks, the differences in WSS driving developmental vascular remodelling. Given the cross-species corroboration of important roles of WSS and RBC dynamics in vascular remodelling through different vertebrate model systems, future work should investigate the relevance of the findings to human physiology and whether the predicted WSS differences are sufficient to trigger pathological vascular remodelling.

Concluding remarks
In summary, our study reports a new mechanism for enhancement of the WSS differences driving vascular remodelling during development. These enhanced differences arise due to the highly heterogeneous distribution of RBCs within the primitive plexus, which is primarily governed by the plasma skimming effect. Additionally, we speculate that vascular remodelling driven by the principle of removing RBC-poor vessels (which are inefficient in transporting oxygen) from the primitive vasculature formed in the course of angiogenesis will ultimately lead to a network layout that avoids portions of the tissue being vascularized but poorly oxygenated. This RBC-driven process, which is highly dynamical and emerging in nature, can importantly contribute to the optimal patterning of vascular networks during development. Conversely, it provides a vascular remodelling mechanism capable of linking changes in RBC royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210113 deformability reported in diseases such as diabetes mellitus or hypertension and the associated microangiopathic complications, which reintroduces the WSS differences in adult microvascular networks akin to those in developmental stage due to abnormal haematocrit distribution arising from locally impeded RBC perfusion. Beyond these findings, our study also has important implications for the mathematical modelling of microvascular haemodynamics in general. In a network of microvessels, multiple effects inexplicable by continuum flow models (albeit widely employed for modelling blood flow in complex vascular networks by existing studies) occur owing to the particulate nature of blood, i.e. essentially a suspension of RBCs. Conventional assumptions such as Poiseuille Law, velocity-/haematocrit-profile symmetry, spatial-/time-average accuracy are all subject to scrutiny when quantifying flow variables such as effective viscosity and WSS in the microcirculation, especially for vessels with a lumen smaller than the undeformed size of an RBC (i.e. roughly 6-8 µm).

Preparation of mouse retina for binary mask acquisition
The mouse strain used in the present study was C57/BL6J. Mice were maintained at the Max Delbrück Center for Molecular Medicine under standard husbandry conditions. Animal procedures were performed in accordance with the animal license X9005/15. Mouse eyes were collected at P5 and fixed with 4% PFA in PBS for 1 h at 4 C, and retinas were then dissected in PBS. Blocking/ permeabilization was performed using Claudio's Blocking Buffer (CBB) [72], consisting of 1% FBS (Gibco), 3% BSA (Sigma-Aldrich), 0.5% Triton X-100 (Sigma-Aldrich), 0.01% sodium deoxycholate (Sigma-Aldrich), and 0.02% sodium azide (Sigma-Aldrich) in PBS at pH 7.4 for 2 h with rocking at 4 C. Primary antibodies were incubated at the desired concentration in 1 : 1 CBB/PBS with rocking at 4 C overnight and secondary antibodies were incubated at the desired concentration in 1 : 1 CBB/PBS for 2 h at room temperature. Retinas were mounted on slides using Vectashield mounting medium (H-1000; Vector Labs).
The Images were taken at room temperature using Zen 2.3 software (Zeiss). Tiled scans of whole retinas were analysed with ImageJ to generate binary masks of ICAM2 and Collagen IV.

Network characterization and region of interest selection
Given the small size of a single RBC relative to the dimension of the whole-plexus vasculature, direct simulation of cellular blood flow in the entire network would be prohibitively expensive if the grid resolution required to resolve the RBC is applied to the flow domain. Therefore, it is necessary to simulate RBC flow in a reduced number of ROIs for investigation of its effect on vessel remodelling.
There are mainly four criteria for the selection of ROIs from the whole-plexus vasculature without losing generality. First, the selected ROIs should be located within the remodelling region of the plexus (which has been estimated to fall into a range of 0.7 times the radial distance between the optic disc and the farthest vessel in the plexus periphery [40]), thus avoiding the sprouting front where vessel anastomoses occurring due to active EC proliferation could present a similar ICAM2 positive and Col.IV negative signature. Second, the selected ROIs should be representative of different zones within the remodelling area, e.g. close to the feeding arteriole, close to the discharging venule, or relatively central in the capillary bed. Third, each ROI should have plenty of vessel regression or stenosis events for targeted observation of the corresponding RBC dynamics. Fourth, each ROI should be of limited size to allow for sufficient steadystate RBC data (i.e. excluding the transient time steps of the simulation until the RBCs have been fully populated with a stable number of cells count from the simulated ROI) with tractable computational cost.

Whole-plexus simulation
A 3D flow model of the luminal surface (electronic supplementary material, figure S1b) is reconstructed from the Col.IV binary mask (electronic supplementary material, figure S1a) using the open-source software PolNet [44], under the assumption of circular vessel cross-sections. The flow domain is then uniformly discretized into cubic lattice grids with a voxel size of Δx. The lattice resolution is determined as Δx = 0.5 µm (corresponding to a simulation time step length of Δt = 4.17 × 10 −8 s) for the mouse retina used here such that the flow can be reliably solved in the reconstructed network (containing 39 514 304 voxels in total). We have previously found that accurate flow rates and WSS values can be recovered if the majority of vessels have a diameter of at least 3Δx and 7Δx, respectively [43].
The whole-plexus simulation in the reconstructed network adopts the non-Newtonian Carreau-Yasuda rheology model (NNCY) [77] following our previous approach [43], where blood is modelled as a homogeneous non-Newtonian fluid with shear-thinning behaviour. By imposing a physiological ocular perfusion pressure (OPP) between the artery and veins of our V-A-V type plexus, a steady flow within the vascular network is solved (electronic supplementary material, figure S1c) to provide boundary conditions for subsequent RBC simulations in designated ROIs of the plexus. OPP = 55 mmHg is chosen for the present simulation based on a literature survey and sensitivity analysis conducted in [43]. Electronic supplementary material, table S1, provides key parameters of the whole-plexus simulation. The computational framework based on immersed-boundary lattice-Boltzmann method is detailed in §S1 of the electronic supplementary material.

RBC simulation in network subsets
To create RBC simulations in ROIs of the retinal network where evident vessel regression events are observed (e.g. electronic supplementary material, figure S1d), we first clip the designated ROI from the whole plexus as a geometric subset (electronic supplementary material, figure S1e). For a ROI subset with N open boundaries, we set up (N − 1) Poiseuille velocity inlets/outlets (where parabolic velocity profiles are imposed with a centreline velocity ofû ) and one pressure outlet (where a reference pressure p out = 0 is set) (see electronic supplementary material, figure S2, for an illustration of the boundary conditions in ROI-1). The velocity boundary conditions are specified based on the NNCY simulation through integrating local flow rates Q at designated vessel cross-sections and subsequently calculating the centreline velocityû from Q under the assumption of Poiseuille flow, followingû ¼ 2 u ¼ 8Q=(pD 2 vessel ). With all the above boundary conditions set up (see electronic supplementary material, table S3), a plasma flow simulation is initiated in the ROI with the fluid viscosity equal to that of plasma η plasma . Once steady flow is achieved and the velocity field is verified against the NNCY simulation, we populate the ROI with RBCs which are continuously fed at a discharge haematocrit of 20% from all inlets of the subset network (electronic supplementary material, figure S1f ). Details of the RBC model can be found in table S2 of the electronic supplementary material.

Vessel selection and quantification of RBC flow
For quantification of RBC perfusion within the selected ROI subsets, we first locate all divergent bifurcations (i.e. composed of one feeding parent branch and two downstream child branches) encountered by the cellular flow in each ROI via exhaustively examining the flow directions in every single vessel segment (ensured identical between the preliminary ROI plasma flow simulation and the whole-plexus simulation) (electronic supplementary material, figure S5a-c). Only vessel branches from these divergent bifurcations are then selected for analysis as representation of the systematic RBC partitioning occurring in the network. In total, 22 divergent bifurcations (excluding bifurcations with excessively short branches or no RBCs involved) and 54 independent vessel segments (excluding 12 repetitive branches) are identified from ROI-2, ROI-3 and ROI-4. To quantify the time-average distribution of RBCs at each divergent bifurcation (electronic supplementary material, figure S5d-f ), we introduce two variables following the practice of [30,49]: Q Ã blood , denoting the proportion of blood flow that a given child branch receives from its parent branch (0 Q Ã blood 1); Q Ã rbc , denoting the proportion of RBC flux likewise (0 Q Ã rbc 1). If Q Ã rbc . Q Ã blood , the child branch is receiving more RBCs than linear allocation and we define it as an 'RBC-enriched' vessel; if Q Ã rbc , Q Ã blood , the child branch is receiving fewer RBCs than the linear hypothesis and we define it as an 'RBC-depleted' vessel.

Evaluation of simulation data against the phase separation model
The PSM proposed by Pries and co-workers [49,51] derived a set of empirical formulations based on experimental observation of arteriolar bifurcations in rat mesentery, and established a flowmediated mechanism to quantitatively describe the RBC fluxes received by child branches of diverging bifurcations within a microvascular network. The PSM correlates the fractional RBC flux FQ E in a child vessel of a divergent bifurcation with the fractional blood flow FQ B that it receives: 1þe À[AþB ln ((FQ B ÀX 0 )=(1À(FQ B þX 0 )))] , X 0 , FQ B , 1 À X 0 0, FQ B X 0 1, FQ B ! 1 À X 0 , 8 > > < > > : (5:1) where A, B and X 0 are fitting parameters derived via linear regression analysis. Physically, A reflects the size difference of the two child vessels, B reflects the shape of the haematocrit profile in the parent vessel and X 0 is related to thickness of the cell-free layer near the corresponding vessel wall. Details of the PSM formulation we have adopted together with the methods and results for evaluating our simulated bifurcations in ROI-2, ROI-3 and ROI-4 against empirical predictions of the PSM are included in §S2 of the electronic supplementary material (including electronic supplementary material, figures S7-S8 and tables S4-S5).

Statistical analysis
Statistical tests are conducted without any removal or modification of the original data. Testing of data normality is performed with the Shapiro-Wilk test and quantile-quantile plots. For normally or approximately normally distributed data, the parametric Welch's T test (without assuming equal variances between two independent samples) is used and the means of data are measured, e.g. figure 3i,j. For asymmetrically distributed data, the non-parametric Mann-Whitney U test is used and the medians of data are measured, e.g. figure 1d,f and figure 2d,e,f,g. All statistical tests are non-paired and two-sided. The difference between groups is considered statistically significant for *p < 0.05 (**p < 0.01, ***p < 0.001, ****p < 0.0001); p > 0.05 is considered non-significant (n.s.).
Ethics. Animal procedures were performed in accordance with the animal license X9005/15. For growing and breeding of transgenic zebrafish lines, we comply with regulations of the ethical commission animal science of MDC Berlin and with FELASA guidelines [74].
Data accessibility. Data and information supporting this article are provided in the electronic supplementary material, including figures S1-S9, tables S1-S5, equations S1-S6 and movies S1-S4. The source code for the blood flow simulations is available at https://github. com/hemelb-codes/hemelb.