Abnormal morphology biases hematocrit distribution in tumor vasculature and contributes to heterogeneity in tissue oxygenation

Significance Oxygen heterogeneity in solid tumors is recognized as a limiting factor for therapeutic efficacy. This heterogeneity arises from the abnormal tumor vascular structure. We investigate the role that anomalies in red blood cell transport plays in establishing oxygen heterogeneity in tumor tissue. We introduce a metric to characterize tumor vasculature (mean vessel length-to-diameter ratio, λ) and demonstrate how it predicts tissue-oxygen heterogeneity. We also report an increase in λ following treatment with the antiangiogenic agent DC101. Together, we propose λ as an effective way of monitoring the action of antiangiogenic agents and a proxy measure of oxygen heterogeneity in tumor tissue. Unraveling the causal relationship between tumor vascular structure and tissue oxygenation will pave the way for new personalized therapeutic approaches.


Supporting Information Text
The Supplementary Information is organised as follows. First, we provide experimental evidence which supports the findings that vessel lengths and diameters are uncorrelated in tumour environments. Next, we describe the fluid structure interaction (FSI) algorithm used for the red blood cell (RBC) simulations and the method used to calculate the width of the cell free layer (CFL). Next, we present our hybrid model of tissue perfusion and introduce our new haematocrtic splitting (HS) model. Finally, we comment on the higher mean oxygen values predicted by our oxygen perfusion model for small λ values.
Vessel lengths and diameters in tumour microvasculature are uncorrelated. In Supplementary Tables S4 and S5 we list Pearson's r-values quantifying the correlation between vessel lengths, L, and diameters, d, [1] where cov(i, j) is the covariance of two variables and σi is the standard deviation of variable i, for the three tumour cell lines used in our experiments. Results are presented for each mouse and each scan. Day 0 was chosen as the day when the tumour vascular network appeared to be fully formed. This typically occurred approximately 8 days after tumour induction, when the tumour size was approximately 4 mm in diameter. We note also that the duration of the observation period is cell-line specific; some tumours grew faster than others and, as a result, soon started pushing on the window, and in such cases the animal had to be culled as per licence limitations. The Pearson's r-values are too low to conclude that a correlation exists between L and d in the tumour vascular networks studied.
Red blood cell suspension model. The lattice Boltzmann method (LBM) numerically approximates the solution of the Navier-Stokes equations for a weakly compressible Newtonian fluid discretised on a regular lattice. We employ the D3Q19 lattice, the Bhatnagar-Gross-Krook collision operator extended with the Guo forcing scheme (1), the Bouzidi-Firdaouss-Lallemand (BFL) implementation of the no-slip boundary condition at the walls (2), and the Ladd implementation of the velocity boundary condition for open boundaries (3). These methods have been extensively used and analysed in the literature (see (4,5) for a detailed description).
The RBC membrane is modelled as a hyperelastic, isotropic and homogeneous material, following the model described in (6). The total membrane energy W is defined by W = W S + W B + W A + W V , where the superscripts denote energy contributions due to strain, bending, area and volume. We employ the surface strain energy density w S proposed by Skalak et al. (7): where κs and κα are the shear and dilation moduli, λ1, λ2 are the local principal in-plane stretch ratios (see (8) for calculation procedure), and W S = dA w S . The shape of the discocyte membrane is approximated by a number N f of flat triangular faces, and W S is numerically calculated based on a finite element method (FEM) approach as where A (0) j is the undeformed area of face j. The bending energy of the RBC membrane is numerically calculated as where κB is the bending modulus, θi,j is the angle between the normals of two neighbouring faces i and j, and θ i,j is the same angle for the undeformed membrane. Finally, we penalise deviations of the total membrane surface area and volume by defining two additional energy contributions: where κA, κV are the surface area and volume moduli, A and A (0) are the current and undeformed membrane surface areas, and similarly with V . The principle of virtual work yields the force acting on each membrane vertex i at position xi through The immersed boundary method (9) is used to couple the fluid and membrane dynamics. The fluid velocity is interpolated at the positions xi of the RBC mesh vertices, and a forward-Euler scheme is used to advect the vertices to satisfy the no-slip condition. The vertex forces Fi are spread to the lattice where they are used as input to the forcing term in the LBM, which ensures local momentum exchange between the membrane and the fluid. See (6) for a detailed numerical analysis of the algorithm.
The RBC model contains five parameters (κs, κα, κB, κA, and κV ). While κs and κB are known from experiments (see review in (10)), the exact values of the three remaining parameters are chosen to ensure that local area, total surface area and volume drift are constrained within a few percent and simulations are stable (see analysis is (6,8)). Table S7 summarises all the parameters in the model.

CFL width calculation.
To calculate the CFL width in channel 1 of the domains in Figure 3, let us consider a vessel cross-section of diameter d at distance l downstream from the first bifurcation in the network. The RBC density, φ(r, θ, l, t), is 1 if there is a RBC at time t occupying the point with radial coordinate r and angular coordinate θ of the cross-section and 0 otherwise. The average RBC density flux Φ(l) going through the cross-section is where v is the fluid velocity, n is the cross-section normal vector and N is the number of simulation time steps in the average (0.5 s of real time simulation sampled every 0.0215 s, N = 23, in our case). We define χ = 0.01 as the average fraction of RBC density flux crossing the CFL. Now we are able to numerically determine the local CFL width W (l, θ): consider a 2D-cone centered and contained in the cross-section with orientation θ and size ∆θ = π/2. The width W (l, θ) is the distance such that Since we are only interested in the spatial evolution of the CFL, the specific value of χ used in the definition is arbitrary. The choice of χ will change the width of the CFL after symmetry recovery, but it will not affect the local characterisation of the CFL spatial evolution after a bifurcation. For example, for any value of χ, the CFL recovery distance can be calculated as the shortest distance l for which the CFL width W do not depend on coordinate θ.

Hybrid model of oxygen transport in vascularised tissue.
Choice of vessel diameters and branching angles in vascular networks. In the branched networks used, we fix the diameter of the inlet vessel so that d inlet = 100 µm. The diameters of the two child vessels (dα and d β ) are assumed to be equal and determined from the diameter of the parent vessel (dP ) via Murray's law (11) so that: Since the network is symmetric about its central axis, vessels on the converging side have the same diameters as those of the same generation on the diverging side (see Supplementary Figure S4a). For all simulations the networks have 6 generations of vessels. The length L of a vessel segment in a given network is related to its diameter d via L = λd, where the positive constant λ is network-specific. For complete specification of the network geometry, in two-dimensional Cartesian geometry, it remains to embed the network in a spatial domain. This is achieved by specifying either the branching angles, or (equivalently) the lengths of the projections of the vessels on the y axis. Denoting by L vert Poiseuille's law and the Fahraeus-Lindquist effect. We simulate flow in the branched networks by following the approach of Pries et al. (12). For blood vessels of length L and diameter d, we assume Poiseuille's law Lµ , [10] where Q is the vessel flow rate, ∆p is the pressure drop along the vessel and µ is the effective viscosity of blood (13). Following (14) we assume that the blood viscosity depends on vessel diameter and haematocrit via the empirical relationship: where µp is the plasma viscosity, H is the vessel discharge haematocrit, Introducing signed flow ratesQi for the sake of brevity, we impose conservation of blood and haematocrit at each network bifurcation, so that iQ i = 0, [11] and iQ iHi = 0. [12] In Eq. (11) and Eq. (12) we sum over the three vessels that meet at that bifurcation. At diverging bifurcations, we impose a HS rule: we use (1) from the main text when CFL memory effects are neglected and (2) from the main text when they are included. Denoting by NB the number of network bifurcations and NV the number of vessels, we have NB unknown pressures P , NV unknown flow rates Q and NV − 1 unknown haematocrit levels (the inlet haematocrit being prescribed)altogether NB + 2NV − 1 unknowns. At the same time, we impose Poiseuille's law (Eq. (10)) for every vessel (NV times), conservation of blood (Eq. (11)) and haematocrit (Eq. (12)) at every bifurcation node (NB times), and an HS rule at all diverging bifurcations (NB/2 times), yielding a total of NV + 5NB/2 algebraic equations. Since every bifurcation connects 3 vessels, we have NV = (3NB + 2)/2, where every vessel is counted twice, except for the inlet and outlet vessels (+2 in the numerator). From this, it follows that the number of equations (NV + 5NB/2) equals the number of unknowns (NB + 2NV − 1).
Oxygen distribution in tissue. In this section, we determine the oxygen concentration c in the tissue. Following (15), we assume that the dominant processes governing its distribution are delivery from the vessel network (via one-way coupling with Eq. (12) and the haematocrit models (1) or (2), i.e. c depends on H l but not vice versa), diffusive transport through the tissue, and consumption by cells in the tissue. We focus on the long time behaviour and, therefore, adopt a quasi-steady state approximation (16)  13) we assume that the oxygen is supplied by vessels to the tissue at a rate which is proportional to their circumference πd l , the vessel permeability γ, where cstp denotes an ideal gas concentration at standard temperature and pressure, p ref denotes the reference partial pressure at the inlet vessel, and α ef f denotes the volumetric oxygen solubility (17). A summary of the parameter values used to solve Eq. (13) is presented in Supplementary Table S8.

Derivation of, and justification for, the HS model with CFL memory.
Parameter dependencies in HS model without memory from (18). The dependencies of the parameters A, B and X0 (see (1) from the main text) on the diameters of the parent and child vessels (dP , dα and d β , respectively), and the discharge haematocrit HP in the parent vessel were first introduced in (12,19) and later adjusted in (18) to achieve a better approximation under extreme combinations of dα, d β , dP and HP . We will use the functional forms from (18), which read B = 1 + 6.98 1 − HP dP [15] and X0 = 0.964(1 − HP )/dP . [16] These functional forms assume that dP is dimensionless and given by dP =d P 1µm , wheredP is the dimensional diameter. We maintain this convention throughout this section.

HS model with memory.
Simplifying assumptions.
Before we explain how we extend the model from (18) to incorporate memory effects, we comment on its main simplifying assumptions. At present, our model does not include any information on local flow rate. Furthermore, the current model does not account for the angle defined by the planes containing the current and previous bifurcation in the network. These simplifying assumptions could easily be relaxed.

Rewriting of the model.
In this section, we rewrite the HS model with memory effects ((2) from the main text) in terms of discharge haematocrit levels H and flow rates Q experienced by the vessels belonging to a given bifurcation ((3) from the main text). The definitions of F Q E,f and F Q B,f can be written as: Substituting these expressions into (2) from the main text gives: Appealing to conservation of blood (overall) QP = Q f + Qu and RBCs (in particular) QP HP = Q f H f + QuHu [17] at diverging bifurcations, we arrive at This equation can also be written as which yields

Choice of parameter values and CFL recovery function.
Here we introduce the functional forms for A f , X 0,f and X0,u, using empirical data to justify our choices. Guided by the dependence of A on the network branching history described in (19) (see Figure 7 therein), we propose where A is given in Eq. (14), the positive constant A shif t corresponds to the maximum CFL disruption effect, and the function f (l; dP ) describes how the recovery of the CFL depends on the distance l to the previous bifurcation and the diameter dP of the parent vessel * . For parameter A, we only have access to the scattered data with respect to the regressor from (12,19) (as opposed to the regressor from (18)), which reads A = −6.96 ln dα d β /dP . [19] Using the extreme values of A in these data (see Supplementary Figure S8c), we estimate A shif t = 0.5. Note that in branching networks with every pair of child vessels having equal radii, both (12) and (18) (18)). * Consistency of the model requires that For simplicity, we model the CFL recovery using an exponential function f (l; dP ) = e − l ωd P , [20] where ω controls the temporal dynamics of CFL recovery. From (20), we note that the CFL width is (approximately) 90% recovered at a distance l90 = 10dP from the previous bifurcation (see also Figure 3g). Accordingly, we choose ω so that Guided by the dependence of X 0,f on flow history described in (19), we propose f (l; dP )) . [21] Assuming, as a first approximation, that X 0,f + X0,u is constant and independent of the distance to the previous bifurcation (see Figure 3g), we define X0,u = X0 (1 + f (l; dP )) . [22] Validation of the HS model with memory. We validate the HS model with memory by comparing its predictions with results from the RBC simulations in the double-t geometry. We assume that all vessels have the same diameter (d = 33 µm), and that the flow rate splits evenly at both bifurcations. If we assume further that the CFL is fully established at the network inlet vessel, H inlet = 20%, then (1) from the main text supplies H1 = H2 = H inlet = 20%. We use conservation of RBCs Eq. (17) and the new HS model (Eq. (3) from the main text) to estimate haematocrit values in the unfavourable and favourable child branches after the second bifurcation (channels 3 and 4, respectively) for varying inter-bifurcation distances δ. The results are summarised in Supplementary Table S9. For δ = 4d, the new HS model predicts haematocrits within 5% relative error of the values calculated from RBC simulations ( Table 2). Given the uncertainty in determining discharge haematocrit in the RBC simulations and given that the new model neglects effects due to asymmetric streamline splitting (21), we conclude that our new model provides a good, leading-order approximation to the effects of CFL disruption on HS. Finally, we compare the CFL evolution dynamics calculated from the RBC simulations (for θ = 0 and θ = π) with those predicted from the proposed evolution of X 0,f and X0,u (Eq. (21) and Eq. (22)). In the absence of a known functional form relating the CFL width W and the minimum flow fraction X0, we define Eq. (23) is based on the diagram in Supplementary Figure S8a [24] We remark that for a well-established CFL (i.e. l → ∞), Eq. (24) Figure S8b). We postulate that this discrepancy is caused by our oversimplification of the relationship between the CFL width and the minimum flow fraction (Eq. (23)). Nevertheless, we can adjust Eq. (24) so that it is consistent with the established CFL width of 1.8 µm by writing [25] In this case the CFL evolution (for θ = 0 and θ = π) follows a trend similar to that observed in our RBC simulations (Supplementary Figure S8b). In particular, our assumption that l90 = 10dP is in good agreement with our simulation results (see dashed line in Supplementary Figure S8b).

Explanation of higher mean oxygen values for small λ.
We observed that CFL disruption effects increase the mean oxygen concentration in the chosen network (Supplementary Figure S6). Here, we provide an explanation of this phenomenon. We define ∆αH = Hα − HP , ∆ β H = H β − HP , [26] where P is the parent branch and α and β are the child branches of any diverging bifurcation. Conservation of blood and RBCs at this bifurcation then yields Qα + Q β = QP [27] Qα (HP + ∆αH) + Q β (HP + ∆ β H) = QP HP . [28] Combining Eq. (28) and Eq. (27) supplies We deduce that, at diverging bifurcations, the haematocrit level in the child branch with higher flow rate deviates less (in absolute value) from the haematocrit in the parent vessel than the branch with lower flow rate. We note further that all paths connecting the inlet and outlet vessels in the direction of blood flow in a given network are topologically and geometrically equivalent. Therefore, heterogeneity in haematocrit splitting arises solely from CFL disruption effects. If haematocrit is elevated in one of the child branches, its impedance will increase, and, as a result, it will receive a lower flow rate.
Combining these two effects, we see that, in the chosen networks, haemoconcentration in any child branch is more significant than haemodilution in its sibling. As a consequence, and given that the strength of the oxygen source term in Eq. (13) is a linear function of H, we observe higher mean oxygen levels when the effects of CFL disruption are taken into account (especially for small λ). Future work will investigate this effect by making source term a function of RBC mass flux (i.e. QH) or relaxing the assumption that the RBCs have infinite oxygen carrying capacity.        Table 1). The missing datum for tumour 3 on Day 3 is due to the laser on the microscope failing during imaging.   Blood flow separation at the two consecutive bifurcations is shown in dotted green, streamlines are sketched with yellow curved arrows, and the CFL recovery on the favourable (unfavourable) side of the parent vessel after the first bifurcation is sketched in red (blue). Whenever F Q B,f < X 0,f (F Q B,u < X0,u), the favourable (unfavourable) branch only draws blood from the CFL and it thus receives pure plasma. (b) Model of CFL recovery as described by Eq. (25) shows similar trends to and is in satisfactory agreement with the CFL width data from RBC simulations in Figure 3g (given the simplifying assumptions). The established CFL width of 1.8 µm chosen by inspection for this particular dataset. (c) Dispersion of values for A (reproduced using Figure 6 from (19)) is used with the regression from (12) to estimate the value of A shif t ≈ 0.5 in Eq. (18), based on deviation from the regression. We assume the CFL disruption to be the primarily cause of this deviation, and thus its maximum (absolute) value should correspond to l = 0 in Eq. (20) (i.e. f = 1 in Eq. (18)).