Vessel compression biases red blood cell partitioning at bifurcations in a haematocrit-dependent manner: implications for tumour blood flow

The tumour microenvironment is abnormal and associated with tumour tissue hypoxia, immunosuppression, and poor response to treatment. One important abnormality present in tumours is vessel compression. Vessel decompression has been shown to increase survival rates in animal models via enhanced and more homogeneous oxygenation. However, our knowledge of the biophysical mechanisms linking tumour decompression to improved tumour oxygenation is limited. In this study, we propose a computational model to investigate the impact of vessel compression on red blood cell (RBC) dynamics in tumour vascular networks. Our results demonstrate that vessel compression can alter RBC partitioning at bifurcations in a haematocrit-dependent and flowrate-independent manner. We identify RBC focussing due to cross-streamline migration as the mechanism responsible and characterise the spatiotemporal recovery dynamics controlling downstream partitioning. Based on this knowledge, we formulate a reduced-order model that will help future research to elucidate how these effects propagate at a whole vascular network level. These findings contribute to the mechanistic understanding of haemodilution in tumour vascular networks and oxygen homogenisation following pharmacological solid tumour decompression.

The tumour microenvironment (TME) is abnormal and associated with tumour tissue hypoxia [1], which is 2 a known biomarker for poor prognosis [2]. In addition, tumour hypoxia is a source of immunosuppression 3 [3] and constitutes a barrier to the success of recent promising immunotherapeutic approaches [3,4]. 4 As such, researchers have proposed normalising the TME to improve oxygenation and overcome these 5 limitations [3]. One of the abnormalities of the TME is vessel compression [5,6] which is a consequence of 6 the proliferation of cells within a solid tumour and the growth of the tumour against its surroundings [7]. 7 Fang et al. showed that the presence of compressed microvessels in tumor tissue is a promising prognosis 8 predictor in non-small cell lung cancer patients [8]. Furthermore, Chauhan et al. demonstrated that 9 normalising the TME through decompressing solid tumours leads to increased survival rate in animal 10 models as well as increased tumour tissue oxygenation and homogeneity [9]. However, our knowledge of 11 the biophysical mechanisms linking tumour decompression to increased tumour oxygenation is limited. 12 Oxygen binds to haemoglobin in red blood cells (RBCs) and is transported through the vasculature 13 with the RBCs. Early investigations demonstrated that haematocrit (the volume fraction of RBCs in 14 total blood) can vary between child branches of a vascular bifurcation as a function of haemodynamic 15 and geometrical vessel properties (see [10] for a review). Our recent work has shown that vascular 16 development works to avoid vessel segments with rare or transient RBC flow through them [11]. However, 17 the abnormal TME is associated with marked heterogeneity in haematocrit [12], including reports of 18 plasma-only channels [13]. We recently identified reduced inter-bifurcation distance as a source of 19 haematocrit variation via its impact on RBC partitioning at bifurcations [14]. However, the impact that 20 other tumour vascular phenotypes such as vessel compression play on this process is not known.
In vitro work has provided some evidence of how RBC suspensions behave in the presence of a 22 compression in single, straight, channels. At low haematocrits (≤5%), the narrowing of a channel leads 23 to a focussing of the RBCs towards the centre of the channel [15,16]. Fujiwara et al. performed a similar 24 experiment in an asymmetric geometry at higher haematocrits (≤20%) and similarly saw a focussing of 25 the RBCs towards the channel centre. They observed that this is more pronounced at lower haematocrits 26 [17]. The focussing of RBCs is identified to be due to an increased shear rate within the compressed 27 section of the channel [15,17]. However, whether the narrowing in RBC cross sectional distribution post 28 compression is permanent or not is unclear, nor is it clear how the narrowing of the RBCs toward the 29 channel centre has an impact on RBC partitioning at a downstream bifurcation. 30 Investigating the transport dynamics of oxygen and other blood-borne solutes in realistic tumour 31 networks is challenging due to the limited experimental tools available. Several groups have proposed 32 the use of mathematical modelling to bridge this gap [18][19][20][21]. In this work, we propose a computational 33 model to study how vessel compression impacts the partitioning of RBCs at a downstream bifurcation. 34 We report the novel finding that, below a critical haematocrit threshold, vessel compression alters the 35 partitioning of RBCs at the downstream bifurcation due to a change in the cross-sectional distribution of 36 the RBCs induced by the compression. Furthermore, we show that this is independent of flow rate and 37 compression asymmetry. In addition, we report, for the first time, the mechanism and length scale for the 38 cross-sectional distribution of RBCs to return to their pre-compression configuration. Finally, we propose 39 a reduced-order model to calculate RBC partitioning at a bifurcation downstream of a compression in a 40 computationally efficient manner. Future investigations can use this reduced-order model to link vessel 41 compression to tumour tissue oxygen heterogeneity on a whole vascular network level. 42 Taken together, our findings suggest that: a) protection against abnormal partitioning at bifurcations 43 due to naturally occurring morphological variations in vessel cross-section can be achieved during 44 development by homogenising and increasing the average haematocrit in networks, b) increased perfusion 45 (in terms of total flow rate through the network) is not sufficient to reverse the anomalous RBC 46 partitioning due to vessel compression, and c) the link between solid stress and oxygen heterogeneity in 47 tumours, and its reported reversal via stress alleviation, can be partially explained via anomalous RBC 48 partitioning at bifurcations due to compressed vessels.

50
Physical model 51 We model blood flow as a suspension of deformable RBC particles in a continuous plasma phase. The 52 plasma is treated as a continuous Newtonian fluid, with the non-Newtonian properties of blood arising 53 from the presence of the deformable RBCs.

54
The model for the RBC membrane is hyperelastic, isotropic and homogeneous with a membrane 55 energy W = W s + W b + W a + W v [22]. W s is based on Skalak's model for the surface strain energy 56 density [23], W b is a bending energy, W a and W v are penalties on changes in membrane surface area 57 and volume, respectively [22]. We determine the deformability of the RBC by the capillary number, where µ is the fluid dynamic viscosity, γ is a characteristic shear rate, r is a characteristic length (the 59 RBC radius) and κ s is the strain modulus of the RBCs. The capillary number is set to 0.1, unless stated 60 otherwise. The bending modulus is determined from a Föppl-von Kármán number of 400 for a healthy 61 RBC [24]. The radius of the RBC is set to r = 4 µm.

62
Numerical model 63 We solve the fluid model with the lattice Boltzmann method (LBM) and the deformable RBC model 64 with the finite element method (FEM). The fluid structure interaction is modelled with the immersed 65 boundary method (IBM). This is a previously validated model from our group [22,25]. The numerical 66  [26], and the 67 simulations have been run on the ARCHER supercomputer.

68
Our LBM algorithm employs the D3Q19 lattice, the Bhatnagar-Gross-Krook collision operator, Guo's 69 forcing scheme [27], the Bouzidi-Firdaouss-Lallemand no-slip boundary condition at the walls [28], and 70 the Ladd velocity boundary condition for inlets/outlets [29]. 71 We discretise the RBC membrane into 720 elements and use a lattice voxel size of 0.6667 µm for 72 the flow simulation, which is sufficient to resolve both the flow and the RBC membrane dynamics [22, 73 30]. We set the viscosity of the plasma and fluid inside the RBC to 1 mPa s. The density of the fluid is 74 1000 kg/m 3 . We chose a dimensionless lattice-Boltzmann relaxation time of 1, which gives a time step 75 of 7.41 × 10 −8 s.

77
We produced four geometries representing a vessel with diameter D and a downstream bifurcation. The 78 first geometry is a control (no compression, Figure 1a), while the remaining three geometries feature a 79 single compression upstream of the bifurcation. The first compression model ( Figure 1b) contains a long 80 compression without a recovery length between compression and bifurcation. In the second compression 81 model (Figure 1c), the compression is short, and there is a short recovery length between compression 82 and bifurcation. The last geometry ( Figure 1d) features a short compression followed by a long recovery 83 segment. The relevant geometrical parameters of all four geometries are summarised in Table 1. 84 We set the channel diameter to D = 33 µm, a typical value for the tumour microvasculature [12]. 85 We assume the cross section of the channel to be circular, except for the section that is compressed, 86 where it takes an elliptical form. We assume the perimeter of the cross section to be constant along 87 the channel, setting the ellipse perimeter to the same value as the uncompressed circular cross section. 88 The segment with elliptical cross section has an aspect ratio of 4.26 [8]. The assumption of an elliptical 89 cross-section within the compression is in line with observations from tumour histological slices where 90 vessel compression is commonly reported as the aspect ratio of the elliptical shape of the vessel cross 91 sections [8,[31][32][33].

92
Our aim is to focus on the effect of the compression. Therefore, we remove any effect from the slope 93 leading to the compression by having a steep transition to and from the compression. We also remove 94 the effect of a bifurcation asymmetry by having both child branches at the same diameter and angle 95 from the parent branch.

96
Inlet and outlet boundary conditions 97 We set the outflow boundary conditions at the child branches as a Poiseuille velocity profile with an 98 imposed maximum velocity and control the ratio of these velocities such that one child branch receives 99 80% of the flow and the other child branch 20%. Unless specified otherwise, the inlet branch has an 100 average velocity of 600 µm/s, a typical value for the tumour microvasculature [12].

101
The inlet boundary condition is an arbitrary pressure value that has no impact on the simulation 102 result. In order to reduce any memory effects and establish a quasi steady-state distribution of RBCs, 103 the cells flow through a straight tube with a length of 25 tube diameters before entering the compression 104 [34]; we call this length the initialisation length. 105 We vary the value of haematocrit within our system from 10% to 30%, covering a wide range which 106 is physiologically present within the tumour microvasculature [12].

107
The physical Reynolds number of the system is 0.04. Therefore, viscous forces dominate and the 108 system is in the Stokes flow regime. For computational tractability we set the numerical Reynolds 109 number to 1, where inertial forces still do not play a significant role.

110
Processing results

111
All haematocrits reported are discharge haematocrits. We calculate the discharge haematocrit, H d , by 112 calculating the fraction of RBC flow to total blood flow at any channel cross section normal to the 113 direction of flow, where Q RBC is the volumetric flow rate of RBCs and Q blood is the volumetric flow rate of blood, 115 i.e. plasma and RBCs. The flow rate of RBCs is calculated by counting the RBCs crossing a plane 116 normal to the direction of flow over a given period of time, ∆t. Knowing the volume of an RBC, 117 V RBC = 100 µm 3 , one can calculate the RBC flow rate, where N is the number of cells that have crossed the plane.

119
In order to quantify the distribution of the RBCs in a cross section, we measure the root mean 120 squared distance (RMSD) of the RBC centres of mass with respect to the channel centreline, where x i is the i th RBC position, and x 0 is the channel centre, both taken on a cross section normal to 122 the direction of flow at points of interest. Figure S1 illustrates how we obtain the positions x i in practice. 123 We non-dimensionalise length by the vessel diameter D, unless stated otherwise. By definition, we 124 set the downstream end of the compression as the reference point with an axial position of 0. Axial 125 positions are positive in downstream direction (l > 0) and negative in upstream direction (l < 0).

126
The separatrix is an imaginary surface separating fluid particles going to one child branch from those 127 going to the other child branch. It is an important tool for the investigation of RBC partitioning at 128 a bifurcation [35]. We determine the separatrix by completing a simulation without RBCs to obtain 129 streamlines unperturbed by RBCs (see Figure S2 for details).

131
At 10% haematocrit, vessel compression alters RBC partitioning at a down-132 stream bifurcation 133 We start by investigating whether vessel compression has an impact on RBC partitioning at a downstream 134 bifurcation. Simulation results for a long compression at 10% haematocrit ( Figure 1b) reveal that the 135 RBC split at the bifurcation is strongly affected by the compression. Figure 1e shows that, for a long 136 compression, the child branch with the lower flow rate is almost depleted of RBCs and has approximately 137 0.5% haematocrit, whereas the control simulation indicates that the same branch has ∼8% haematocrit 138 in the absence of a compression. The control simulation is in agreement with the standard plasma 139 skimming model, see Figure S3. 140 The short compression short recovery geometry ( Figure 1c) shows a smaller impact on the RBC split 141 than the long compression. With ∼ 3.5% haematocrit in the child branch with the lower flow rate, this 142 is still less than half of the haematocrit in the control simulation. In order to test whether the symmetry of the compression has an effect on the partitioning of the 144 RBCs at the downstream bifurcation, we investigated an asymmetric compression ( Figure S4). The 145 simulation results are similar to those obtained with a symmetric compression. Thus, we infer that, 146 under the present conditions, the asymmetry of the compression does not have an important effect on 147 the downstream partitioning of RBCs. 148 We also investigated the effect of flow rate on RBC partitioning. By changing the flow rate, we change 149 the capillary number, which quantifies the deformation of the RBCs. We performed two additional 150 simulations at a capillary number of 0.02 and 0.5, to cover the RBC tumbling and tank-treading regimes 151 and the range of flow rates typical for the tumour microvasculature [12]. Figure S5 shows that the 152 RBC partitioning does not change with capillary number, which implies that the flow rate and capillary 153 number are not important parameters for RBC partitioning in the presence of a compression within the 154 studied range of flow rates.

155
Narrowing of cell distribution alters partitioning of RBCs 156 Next we investigated which mechanism leads to the observed changes in partitioning, and why the 157 different geometries have different effects. 158

5/21
As blood flows through the compression, the shear rate increases since 1) the fluid velocity within the 159 compression is larger due to mass conservation and 2) the width of the channel along the compression 160 axis is reduced. Our simulations show that RBCs situated close to the wall prior to the compression 161 migrate across streamlines towards the channel centre. After leaving the compression, RBCs do not 162 immediately migrate back towards the wall. As a consequence, the RBC distribution downstream of the 163 compression is more narrow than upstream of the compression. This explanation is in line with prior 164 findings from experimentalists [17]. Figure 2a-d illustrates this mechanism. 165 In order to quantify the narrowing of the RBC distribution, we plot the RMSD of the RBC centres 166 of mass along the compression axis. Figure 2e-g shows that, for all geometries, there is a narrowing of 167 the distribution of cells in and after the compression compared to the region before the compression. 168 However, in the short recovery geometries (Figure 2f-g), the RBC distribution partially recovers before 169 the cells reach the bifurcation. This explains why the RBC partitioning is more affected when there is 170 no recovery length between compression and bifurcation ( Figure 2e). Previous studies report an increase 171 of the cell free layer (CFL) post compression compared to the CFL thickness pre compression, here seen 172 as a narrowing of the RBC distribution, which leads to a partitioning bias of RBCs towards the higher 173 flowing branch [36]. 174 In order to investigate the behaviour of the RBC distribution after the compression, we increased the 175 distance between the compression and the bifurcation from 2D to 25D (Figure 1d). Figure 2g shows a 176 gradual recovery of the RBC distribution between the compression and bifurcation, although 25 channel 177 diameters are not sufficient to reach the same RMSD as before the compression.

178
While our data imply that a mechanism exists that leads to the recovery of the RBC distribution, it is 179 not clear a priori what the underlying mechanism is. We assume that, given enough channel length after 180 the compression, the RBC distribution will fully recover eventually. Katanov et al. demonstrated that, 181 from an initially uniform distribution of RBCs in a channel, the formation of a stable CFL is governed 182 by the shear rate time scale and takes a length of about 25 vessel diameters to form, independently of 183 flow rate, haematocrit or vessel diameter [34]. Our data suggest that the opposite effect, the recovery of 184 an initially heterogeneous RBC distribution where most of the RBCs are close to the channel centre, 185 cannot be described in the same way since a length of 25 channel diameters is not sufficient for recovery. 186 The shear rate is lower and cells move faster near the channel centreline, which should lead to a weaker 187 shear-induced recovery of the cell distribution along a distance of 25D. We hypothesise that cell-cell 188 interactions are the dominant driver for the recovery.

189
Increasing haematocrit reduces bias in RBC partitioning 190 To test the hypothesis that cell-cell interactions drive the recovery of the distorted RBC distribution, we 191 investigate blood flow at an increased haematocrit of 20%. Figure 3 shows that the long compression 192 no recovery geometry still leads to a deviation from the control simulation. However, the deviation is 193 smaller than at 10% haematocrit ( Figure 1). This observation can be explained by Figure 2e-g which 194 reveals that the narrowing of the RBC distribution is less pronounced at higher haematocrit compared to 195 lower haematocrit. We conclude that, as the haematocrit increases, hydrodynamic cell-cell interactions 196 become more relevant, leading to a smaller narrowing of the cell distribution by the compression as well 197 as a faster decay of the narrowed RBC distribution after the compression.

198
For the short compression long recovery geometry, the deviation from the control is almost non-199 existent at H d = 20% (Figure 3). Whilst the control simulation shows 15.6% haematocrit in the lower 200 flowing child branch, the compression merely alters that value to 13.7%. As can be seen from Figure 201 2f, not only is the narrowing of the RBC distribution in the compression smaller, but after exiting the 202 compression the RBC distribution is much closer to its pre-compression counterpart. 203 We also investigated the role of the distance between the short compression and the bifurcation by 204 increasing it from 2D to 25D. Figure 2g shows that the RBC distribution with H d = 20% eventually 205 recovers and goes back to its pre-compression level, contrary to the simulation at 10% haematocrit.

206
The decreasing effect of the constriction on the RBC distribution at increasing haematocrit raises the 207 question whether there is a critical haematocrit above which the RBC partitioning is not modified by the 208 presence of an upstream constriction. To that end, we increased the haematocrit in the parent branch to 209  30% and revisited the long compression geometry that has no recovery length between compression and 210 bifurcation. Figure 4 shows that there is no significant difference in RBC split when compared to the 211 control without compression. We conclude that the critical haematocrit value lies near 30%.

213
A challenge in the theoretical study of RBC transport in networks is computational expense [25, 37, 214 38]. For this reason, several authors have proposed the use of reduced-order models to quantify the 215 partitioning of RBCs at bifurcations, which is key for tissue oxygenation modelling due to the RBC's 216 role as oxygen carrier. The most common model existing for partitioning of RBCs is that presented by 217 Pries et al. [39,40], although others exist, e.g. [41]. 218 We have demonstrated that vessel compression indeed has an impact on the partitioning of RBCs 219 at a downstream bifurcation and that this is not captured by a state-of-the-art reduced-order model. 220 In order to investigate how this effect propagates on a network level, we propose a novel reduced-order 221 model that captures this phenomenon. We make four main assumptions for the reduced-order model: 222 1. An RBC's centre of mass will go to the same child branch as its underlying streamline. Similarly 223 to reports of RBCs crossing the separatrix prior to a bifurcation [35,36,42], in our simulations we 224 observed < 5% of RBCs near the separatrix crossing streamlines. Due to this small fraction, we 225 deem the assumption appropriate. 226 2. A curved separatrix independent of the diameter ratio between the child and parent branch and 227 independent of Reynolds number is used, whereas the curvature of the separatrix generally depends 228 on both parameters [35,43]. Since blood flow in the microvasculature is in the low-Reynolds regime, 229 a small change in Reynolds number has a negligible effect on the separatrix. The diameter ratio 230 has been shown to have a more significant impact, even at a Reynolds number of 0 [43]. However, 231 the impact on the curvature of the separatrix is higher towards the vessel walls, where there are 232 fewer or no RBCs due to the existence of the CFL. Considering the difficulty of parameterising a 233 curved separatrix as a function of diameter ratios, a curved separatrix for a diameter ratio of 1 is 234 used at the expense of a small but acceptable modelling error. 3. The cross-sectional distribution of RBC centres of mass can be approximated by a step function, 236 whereas the distribution profile of RBCs tends to quickly, but not instantaneously, reduce at the 237 edge of the RBC distribution as is often reported [44,45] and indeed observed in our simulations. 238 The step function, however, is a good fit and simplifies the model considerably. We assume the 239 step function to take the shape of an ellipse on any given cross section of the channel (Figure 5a-c).  4. The cross section that determines which child branch an RBC enters is located about 2D/3 243 upstream of the bifurcation. Our simulations show that this is the upstream perturbation length 244 after which the streamlines start to curve in order to enter the child branches. The length 2D/3 245 is similar to that found in other studies [35]. Therefore, the reduced-order model needs to be 246 able to predict the cross-sectional distribution of the RBCs up to a point 2D/3 upstream of the 247 bifurcation.

248
The step function that approximates the RBC distribution in a channel cross section has the form of 249 an ellipse. Therefore, the major and minor semi-axes of this ellipse, a and b, need to be defined. We 250  Table 2 for 20% haematocrit compared to simulation data.
found that the best results are obtained when 251 1. the aspect ratio of the ellipse is determined by the ratio of the RMSD along the width and height 252 directions (where the height direction is the axis of the compression), 253 2. the ellipse encloses 90% of the RBCs' centres of mass.

254
Next, we propose a function that describes the development of the radius ρ, along the compression 255 axis, of the step function along the channel length l between the end of the compression (defined as 256 l = 0) and the point 2D/3 upstream of the bifurcation. The area of the ellipse is given by A = πab 257 where a and b are the major and minor semi-axis of the ellipse for the step function. Once the minor 258 axis b, which is ρ(l), and the aspect ratio (l) = a(l)/b(l) are known, the model can predict the number 259 of RBCs entering either child branch. 260 Our simulations show that there are three key mechanisms governing the lateral RBC distribution 261 when entering and leaving the compression. The first is that the RBC distribution is suddenly narrowed 262 by the compression. Secondly, upon exiting the compression, the RBC distribution sees a quick but only 263 partial lateral recovery due to the expansion of the streamlines. Lastly, cell-cell interactions lead to a 264 slow recovery of the RBC distribution to its pre-compression distribution via cross-streamline migration 265 if given sufficient length. We model the flow expansion of the RBC distribution with a logistic term and 266  Table 3. Absolute differences of the discharge haematocrit in both child branches between the results of the HemeLB simulations and the predictions of the reduced-order model. In each box, the top left is the difference for the higher flowing child branch, whereas the bottom right is the difference for the lower flowing child branch. See Figure 1a- the cross-streamline recovery with an exponential decay: For l = 0, at the downstream end of the compression, the equation returns the ellipse radius inside the 268 compression, which is within 5% of ρ c , due to the second term becoming very small when l = 0, but not 269 vanishing. ρ sl is the change in RBC distribution due to the flow expansion which occurs over a length 270 scale 2l g . The length l s determines the steepness of the slope. ρ d is the change in RBC distribution due 271 to cell-cell interaction, which occurs over a longer length scale l r l g . For long distances, l → ∞, we 272 have ρ(l) = ρ c + ρ sl + ρ d which is the width of the fully recovered RBC distribution and the width of the 273 unperturbed RBC distribution before the compression. If l is not sufficiently large, the RBC distribution 274 is still affected by the compression, and the RBC split in the bifurcation tends to be biased accordingly. 275 Although it may be possible to construct a reduced-order model with fewer parameters, the advantage 276 of Eq. (5) is that all parameters have a physical meaning and can potentially be predicted by separate 277 models. 278 We obtained the numerical values of the parameters in the reduced-order model by fitting ρ(l) to 279 the simulation data for 10% and 20% haematocrit, respectively. The parameters are listed in Table 2, 280 and Figure 5d shows an excellent agreement between ρ(l) and the simulation data at H d = 20%. In 281 particular, we find that l g ≈ 0.66D, independently of the chosen haematocrit. This value corresponds to 282 the characteristic length which describes the streamline recovery after a distortion. We find that the 283 cross-streamline recovery length l r reduces by a factor of about 2 when the haematocrit is increased from 284 10% to 20%. This is in line with literature reporting that shear-induced diffusion is directly proportional 285 to the particle concentration [46].

286
With the reduced-order model being calibrated, we can now predict the RBC partitioning at the 287 downstream bifurcation and compare these results with actual RBC simulation data. We apply the 288 separatrix model to the cross-sectional RBC distribution predicted by the reduced-order model at the 289 length l that marks the distance between the compression and the point of bifurcation. 290 Table 3 compares the absolute difference in discharge haematocrit obtained from the HemeLB 291 simulations and the reduced-order model. We find that the reduced-order model accurately predicts the 292 impact of the compression on the RBC partitioning at the downstream bifurcation within 1% on average. 293 Notwithstanding the assumptions underlying the reduced-order model, the relative error is low. Our 294 novel approach, therefore, provides a means of modelling the disturbance caused by a compressed vessel 295 in network simulations, which has not been possible using established empirical models [39,41]. Despite 296 this success, further simulations are necessary to extend the validity of the reduced-order model to a 297 larger parameter space. At 10% haematocrit, a converged suspension of RBCs requires a long devel-299 opment length 300 We observed that, at 20% haematocrit, the RMSD of the RBCs after 25D downstream of the compression 301 has recovered to 98% of its original value prior to the compression (Figure 6a). However, at 10% 302 haematocrit in the same geometry, the RMSD recovery is incomplete after the 25D. In fact, the 303 reduced-order model predicts a partial recovery of the RMSD at 50D to only 91% of its pre-compression 304 value (Figure 6b).

305
Given the results in Figure 6b, we hypothesise that, at 10% haematocrit, details of the RBC 306 initialisation in the simulation play a role. An inconsistent RBC distribution upstream of the compression 307 might affect the overall outcome of the simulation. To confirm this, we increased the length of the 308 periodic tube that is used to generate the RBC distribution fed into the compression geometry from 309 25D (Figure 6b) to 100D (Figure 6c). We observed that the longer tube leads to a narrowing in the 310 pre-compression distribution of the RBCs (Figure 6c). At 10% haematocrit, when the initialisation 311 length is 100D, we can assume that the RBC distribution has reached a steady state. In fact, Figure 6c 312 shows that after 50D the RMSD of the RBCs has recovered to 98% of its pre-compression value. Despite 313 the sensitivity of the pre-compression distribution on the cell initialisation strategy at H d = 10%, we 314 found that the RBC dynamics after the compression is quantitatively and qualitatively similar for both 315 RBC initialisation lengths used.

317
The tumour microvasculature is abnormal and linked to tumour tissue hypoxia [1], which is a known 318 biomarker for poor prognosis [2] and a barrier to recent promising immunotherapeutic approaches [3]. 319 One such abnormality is vessel compression [5]. Previous studies have shown that decompressing tumour 320 vessels leads to increased survival rates [8] via increased perfusion [1,8] and oxygen homogenisation 321 [9]. However, the mechanism linking tumour decompression to increased oxygen homogeneity is unclear. 322 Oxygen binds to haemoglobin in red blood cells (RBCs) and is transported through the vasculature with 323 the RBCs. We recently identified the reduced inter-bifurcation distance associated with the pro-angiogenic 324 tumour environment as a source of oxygen heterogeneity via its impact on RBC splitting at bifurcations 325 [14]. However, the impact that other tumour vascular phenotypes, such as vessel compression, play on 326 this process is not known.

327
Motivated by the limitations on experimental methods available to query this process, we propose a 328 computational model to elucidate the link between vessel compression and abnormal RBC partitioning 329 12/21 at bifurcations. Our numerical simulations show that a vessel compression enhances the disproportional 330 partitioning of RBCs at a downstream bifurcation in favour of the higher flow rate child branch. This 331 is a consequence of the previously identified narrowing of the RBC distribution within the vessel cross 332 section [17]. 333 Similarly to previous studies [15,17], we identify the mechanism leading to this narrowing as RBC 334 cross-streamline migration towards the vessel centre due to an increased shear rate within the compression. 335 Once the RBCs leave the compression, their cross-sectional distribution gradually goes back to their 336 pre-compression configuration in a haematocrit-dependent manner. This process is significantly slower 337 at 10% haematocrit, where the dynamics occur over a length of ∼ 50 vessel diameters. However, at 338 20% haematocrit, there is an almost instantaneous recovery, and after 2 vessel diameters in length, the 339 difference in RBC partitioning compared to a control simulation is negligible. Furthermore, we show 340 that at 30% haematocrit, the difference between the compressed and control simulation is negligible. 341 This suggests the presence of a critical haematocrit above which vessel compression no longer alters the 342 partitioning of RBCs at a bifurcation. We hypothesise that the different dynamics at 10%, 20% and 30% 343 haematocrit are caused by cell-cell interaction which increases with haematocrit. We also show that an 344 asymmetric compression does not lead to a measurable difference in the partitioning of RBCs compared 345 to a symmetric compression. Likewise, a reduction and increase in flow rate by a factor of 5, respectively, 346 does not significantly change the RBC partitioning. 347 We propose a reduced-order model to approximate the partitioning of RBCs at a bifurcation 348 downstream of a compression, which we show has an error of 1% compared to our fully resolved 349 numerical simulations, in the range of parameters studied. This model has the potential to overcome 350 the computational tractability limitations associated with simulating RBC flow in large computational 351 domains and will contribute to unravelling the dynamics of oxygen transport in large vascular tumour 352 networks.

353
The implications of our findings are multiple. First, our results show that the effect of vessel 354 compression on the downstream partitioning of RBCs is only apparent when discharge haematocrit and 355 the distance between the compression and the bifurcation is sufficiently low. Another study by our group 356 showed that two consecutive bifurcations within a short distance can also alter the partitioning of RBCs 357 at the downstream bifurcation [14]. Furthermore, we showed that interbifurcation distances are much 358 reduced in the tumour micro-environment, and Kamoun et al. showed that haemodiluted vessels are 359 more common and are present across a larger range of diameters in tumour networks than in controls 360 [12]. Taken together, this suggests that healthy vascular networks are structurally adapted to protect 361 themselves from mechanisms leading to RBC transport heterogeneity and that this may be compromised 362 in diseased networks.

363
Second, previous studies have shown that decompressing tumour vessels leads to increased survival 364 rates [8,9]. This has been attributed to a) increased tumour perfusion due to reduced vessel resistance [1, 365 8] and b) reduced hypoxia fraction and increased oxygen homogeneity [9]. Our results of anomalous RBC 366 partitioning being unaffected by increases in flow rate support the view that increasing total perfusion 367 through the network may not be sufficient to homogenise oxygenation if it is not accompanied by vessel 368 decompression (or other forms of structural remodelling normalising RBC partitioning).

369
Third, Kamoun et al. found that up to 29% of tumour vessels in an animal model of glioma 370 experience haemodilution (defined as having haematocrits below 5%) and proposed a mechanism whereby 371 extravasated plasma from leaky vessels would be reabsorbed by other vessels and lead to haemodilution 372 [12]. Along similar lines, recent studies have reported findings of tissue hypoxia near perfused vessels 373 [2]. Our results demonstrate that, in the presence of vessel compression and uneven flow split at 374 bifurcations, haematocrit can decrease from 20% to nearly 0% following two consecutive bifurcations 375 without contributions from interstitial fluid. Our present findings, therefore, provide an alternative 376 explanation of the occurrence of haemodilution in tumour networks. Future work should elucidate the 377 relative importance of these two mechanisms.

378
Lastly, we identified that, in the semi-dilute regime of 10% haematocrit, achieving convergence in an 379 RBC suspension that has been disturbed requires longer distances than previously thought. Katanov et 380 al. reported that it takes a length of 25 diameters for the CFL of a randomly initialised suspension in a 381 straight channel to converge [34]. However, our data suggest that up to 100D of length is required for an 382 initially compacted RBC distribution to expand and reach a steady distribution when the haematocrit 383 is 10%. This finding supports the view that the cross-sectional distribution of RBCs at low in vivo 384 haematocrits may be away from equilibrium not only in diseased vascular networks, as we previously 385 showed in tumours [14], but also under physiological conditions where inter-bifurcation lengths average 386 fewer than 100D. Further research into the network-level dynamics arising from our results is warranted. 387 As a final comment, our results should be considered for in vitro experiments that need to carefully 388 consider the design of microfluidic devices if full convergence of RBC suspensions is required in the 389 semi-dilute regime. It is also an additional challenge for in silico studies with open boundary conditions, 390 where not only the insertion of cells needs to be considered [47], but also their cross-sectional distribution. 391

392
In this work we have demonstrated that vessel compression can alter RBC partitioning at a downstream 393 bifurcation. Interestingly, this happens in a haematocrit-dependent and flow rate-independent manner. 394 We argue that these findings contribute to the mechanistic understanding of haemodilution in tumour 395 vascular networks and oxygen homogenisation following pharmacological solid tumour decompression. 396 Furthermore, we have formulated a reduced-order model that will help future research elucidate how 397 these effects propagate at a whole network level. Unravelling the causal relationship between tumour 398 vascular structure and tissue oxygenation will pave the way for the development of new therapeutic 399 strategies. Figure S1. (a) Three RBCs before and after a plane of interest. Lines indicate RBC trajectories, assumed as straight. (b) Side view as each cell crosses the plane at a given coordinate (x, y, z). The RMSD is calculated in the compression axis (here seen as height of channel) by setting x 0 as the channel centreline (always zero) and x i as the height coordinate of the RBC as it crosses the plane. This measures the distribution in the height of the channel. For illustration purposes only three cells are shown, whereas several hundred are accounted for. Figure S2. Intuition for separatrix. Blue/red lines are streamlines ending in the top/bottom child branch, respectively. The separatrix is the surface separating the blue from the red streamlines on the plane.