# Spatio-temporal instabilities of blood flow in a model capillary network

Mathieu Alonzo, Nathaniel Karst, Thomas Podgorski, John Geddes, Gwennou Coupier

► To cite this version:

Mathieu Alonzo, Nathaniel Karst, Thomas Podgorski, John Geddes, Gwennou Coupier. Spatio-temporal instabilities of blood flow in a model capillary network. Physical Review Fluids, 2024, 9 (10), pp.104401. 10.1103/PhysRevFluids.9.104401. [hal-04749016]

HAL Id: hal-04749016

https://hal.science/hal-04749016v1

Submitted on 22 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Spatio-Temporal Instabilities of Blood Flow in a Model Capillary Network

Mathieu Alonzo,¹,² Nathaniel J. Karst,³ Thomas Podgorski,²,\* John B. Geddes,⁴ and Gwennou Coupier¹,†

¹Université Grenoble Alpes, CNRS, LIPhy, F-38000 Grenoble, France

²Université Grenoble Alpes, CNRS, Grenoble INP, LRP, F-38000 Grenoble, France

³Babson College, Wellesley MA, USA

⁴Olin College of Engineering, Needham MA, USA

(Dated: October 8, 2024)

We present experimental evidence of multiple blood flow configurations in a relatively simple microfluidic network under constant inlet conditions. We provide evidence of multistability and unsteady dynamics and find good agreement with a theoretical one-dimensional advection model for blood flow in microvascular networks that relies on the widely used laws for rheology and phase separation. We discuss the ramifications for microfluidic experiments and measurements using blood and implications for *in vivo* microcirculation. Our findings suggest that further modeling in microvascular networks should discard the usual assumption of unique, steady-state flow solutions, with crucial consequences regarding gas, nutrient, and waste transport.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Microfluidic system, 2-bridge model network and notations. (a) The channel inlets (bottom-left corner) allow control of the feed haematocrit $H_D$ by mixing of a buffer solution and a concentrated RBC suspension, as well as the inlet average velocity $V_{in}$ that was determined in the depicted zone of interest. For large width or low haematocrit (here $w = 36 \ \mu m$, $H_D = 0.2$), the flow splits symmetrically at the first bifurcation and no flow is observed in the bridges ($V_a = V_b = 0$), similar to an homogeneous Newtonian fluid, while in (b), for thin enough channels and large haematocrit (here $w = 17 \ \mu m$, $H_D = 0.55$), different flow patterns emerge: a symmetry-breaking takes place with non-zero RBC velocities and concentrations in the bridges.</figcaption></figure>

INTRODUCTION

Multiphase fluid systems involving non-trivial rheology are often a source of nonlinear phenomena in network flows [1] as is also the case in models of urban traffic even with idealized configurations [2]. At the microscale, the flow of droplets through networks exhibits bistability and oscillations [3–6]. In the field of biological fluids, oscillations of microvascular flows have been observed *in vivo* [7, 8]. Although active regulation mechanisms such as vasomotricity [9] are involved in local fluctuations of blood flow, the intrinsic stability of a passive microvascular network is an interesting question to address as it might be related to flow configuration in a regulated network.

The human blood network architecture is characterized by the superposition of two components: a tree-like network of arterioles and venules that connects to a dense, mesh-like, capillary bed with diameters ranging from 5 to 100 µm, i.e. roughly 1-10 cell diameters [10, 11]. The multiplicity of paths of various lengths across an organ leads to strong spatial heterogeneities in red blood cell (RBC) and fluid travel times. These various path lengths might lead to heterogeneous oxygenation, and possibly accumulation of various substances involved in certain diseases [12].

Blood is a dense suspension of cells, mainly RBCs and a small fraction of white cells and platelets, suspended in plasma. Because of the complex mechanical properties of blood cells and their interactions, blood dynamics at the cellular scale is inherently fluctuating, either due to RBC shape dynamics [13–15], hydrodynamic interactions and collisions [16] or aggregation [17, 18]. While these small scale fluctuations may eventually induce perturbations of the flow distribution at the network scale, it is of fundamental interest to first consider one-dimensional (1D) models where properties (haematocrit, velocity) have been averaged over the individual vessel’s cross section as an initial approach to blood traffic in capillary networks. In this approach, network traffic is modeled by coupling pressure-flow relationships in network branches, which can be described by e.g., empirical laws for the effective viscosity as a function of local haematocrit [19] and mass conservation and phase separation laws at bifurcations [20–22] that lead to heterogeneity of the haematocrit distribution [23, 24].

The stability of such a strongly coupled network flow problem has been theoretically addressed in the literature [25–32], revealing multistability and sustained oscillations. Despite these predictions, most recent works on microvascular flows still implicitly postulate the existence of a single steady state [12, 33–37], without considering the experimental evidence of possible perturbations from the presence of (rare but larger) white blood cells [38], local adhesion or aggregation [39], and pathological or impaired cells [40]. Notably, no experimental study demonstrating the intrinsic multistability and nonlinear dynamics of blood flow in networks under controlled and steady conditions is reported in the literature.

In this work, we present experimental evidence of the emergence of multiple flow patterns in a network and show that their main characteristics are in agreement with a theoretical analysis based on a 1D advection model, without needing to invoke cell-scale fluctuations. To that aim, we analyse the local evolution of velocity and haematocrit distributions in a simple 2-bridge ladder-like symmetric network (Fig. 1), exploring a wide range of inlet haematocrit and various channel widths.
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Experimental velocity ratios $V_2/V_1$ and $V_3/V_2$ as a function of normalized time, for three channel widths. The inlet haematocrit is constant in time ($H_D \simeq 0.65 \pm 1\%$, see Fig. 9(a) in appendix ) and is chosen large enough to highlight clear behaviours out of the measurement noise. Concentrations in each bridge, normalized by the maximum concentration, are shown in gray scale. Colored bands identify the reference state, accounting for measurement fluctuations. Other acquisitions for the $w = 11~\mu m$ network are shown in Fig. 9(b) and exhibit similar non-trivial dynamics. (b) The same data complemented by two other time series for $w = 11~\mu m$ and one time series for $w = 17~\mu m$, plotted in the velocity ratio space. “Up” or “Down” couples in each quadrant refer to the direction of flow in the bridges.</figcaption></figure>

## EXPERIMENTAL METHOD

We detail here the experimental methods, leaving the details to the appendix. The microfluidic chip was fabricated using classic soft-lithography techniques by casting polydimethylsiloxane on a SU-8 master produced by direct laser lithography (Kloé, Dilase 250), and then bonded to a plasma-treated glass slide. All channels in a given 2-bridge network have equal length $L = 0.5$ mm, equal height $h = 30~\mu m$, and equal width $w \in [11, 36]~\mu m$. Blood samples were collected by *Etablissement Français du Sang* in citrate tubes to prevent coagulation. RBCs were separated by centrifugation and re-suspended in a buffer solution made of 68.5% of nominal phosphate-buffered saline (PBS) mixed with 31.5% Optiprep (v/v) and 2 g/L Bovine Serum Albumine. This follows a recently standardized protocol that ensures RBC conservation while also preventing sedimentation [41].

The inlet haematocrit $H_D$ was continuously adjusted by varying the relative flow rates of two separate inlets, achieving a mixing of concentrated RBCs with discharge (reservoir) haematocrit $H_{D0} \simeq 60-70\%$ and the buffer solution (Fig. 1(a)). A critical aspect of the experimental procedure is the ability to quickly rinse the system with buffer solution immediately after an acquisition, at constant wall shear stress, thus confirming the absence of clogging and validating that observed phenomena are intrinsic to the nonlinear dynamics arising from rheology/phase separation coupling (see Movie in [42]).

The discharge haematocrit $H_D$ in the inlet branch, i.e., the ratio of RBC flow to total flow, was determined by measuring the tube haematocrit by light absorption [16, 43], which is then corrected for the Fåhræus effect [44]. A PIV-like method based on the intensity patterns of the RBCs was developed to measure the maximum (central) velocity in the inlet ($V_{in}$, of the order of the mm/s) and in the main branches ($V_1$ to $V_3$, see Fig. 1). In the following, time is rescaled by the bulk transit time $t_0 = 8L/\langle V \rangle$, where $\langle V \rangle$ is the mean inlet velocity, approximated to $V_{in}/2$ in the experiments. Deviations of the velocity ratios $V_2/V_1$ and $V_3/V_2$ from 1 indicate the existence of flux in the bridges; this may happen even if no cells enter the bridge, since a low flow in the bridge will be fed only by the cell-depleted fluid flowing near the wall of the upstream branch [45, 46].
MULTIPLICITY OF FLOW SOLUTIONS

The flow solution for a homogeneous, Newtonian fluid flowing in our model symmetric network, hereafter called the reference state, is trivial with $V_a = V_b = 0$. Here, when channels are thin enough and haematocrit is above a threshold, the observation of the network over long times with fixed inlet conditions reveals a rich dynamics featuring successions of apparently stable asymmetric states and quick spontaneous transitions involving flow reversal in the bridges (Fig. 2(a)). Noteworthy, some states may persist over extremely long times compared to the typical transit time $t_0$. For the thinner channel, Fig. 2(b) shows that velocity ratios fluctuate between 0.8 and 1.2, indicating non-negligible flux in both bridges. Because of slight imperfections in the network resulting from the manufacturing process, the reference state does not correspond to velocity ratios strictly equal to 1. Once identified (see appendix ) this reference state is used to normalize results for comparison between different networks.

Fig. 3(a) reports such a comparison while also considering the influence of the feed haematocrit $H_D$. It displays a plot of the velocity ratio $V_2/V_1$ as a function of $H_D$ for a single acquisition beginning about 30 s ($\sim 10t_0$) after blood injection in the network. For narrow enough channels ($w \le 20 \mu m$), Fig. 3(a-i) features two regions: a plateau for $H_D \lesssim 0.45$ (tube haematocrit $H_T \lesssim 0.33$ for $w = 11 \mu m$) followed by an apparently multi-valued scatter, which is absent for larger channels (Fig. 3(a-ii)). These results demonstrate that the existence of multiple nontrivial states is triggered by crossing thresholds in channel width and inlet haematocrit.

Due to the fixed measurement time, these states can be either transient or converged. Because of the long residence times in different states, an extensive and statistically complete study of the nature of these states would require long-lasting experiments (of the order of an hour) that would most probably be interrupted by unavoidable clogging events in such *in vitro* conditions. This experimental constraint makes a thorough experimental study beyond the scope of this article. We show below that the experimental evidence of a rich dynamics characterized by transitions between multiple states is well supported by theory.

MATHEMATICAL MODEL

We base our theoretical analysis on a previously established one-dimensional advection model for the transport of haematocrit through the network [47]. The complex nature of blood flow reveals this model to be fundamentally nonlinear: red cell separation at diverging junctions depends on the flow, diameter, and haematocrit of surrounding vessels; and blood viscosity depends on red cell concentration. Resistance to flow, which determines the advection velocity in every vessel, therefore depends on the haematocrit distribution throughout the entire network. As a result, there may exist multiple equilibrium solutions for flow and haematocrit in the network, and these solutions may themselves lose stability to periodic oscillations [48].

More precisely, if we normalize by the inlet flow, the network under consideration has 9 distinct flows $Q_i$, which are constrained by conservation of mass at each of the 6 nodes in the network. We close the system by identifying three Kirchoff equations $\sum_i Q_i R_i = 0$, where $R_i$ is the resistance in vessel $i$ and is given by Poiseuille's law,

$$R_i = \frac{128L}{\pi w^4} \eta_i, \qquad (1)$$

where $L$ and $w$ are the length and diameter of the vessel, respectively, and $\eta_i$ is the viscosity. We use a well-known empirical model for *in vitro* viscosity as a function of both haematocrit and the vessel's geometric properties [19, 22]:

$$\eta_{\text{vitro}} = 1 + (\eta^* - 1) \frac{(1 - H_D)^C - 1}{(1 - 0.45)^C - 1} \qquad (2)$$

$$\eta^* = 220e^{-1.3w} + 3.2 - 2.44e^{-0.06w^{0.645}} \qquad (3)$$

$$C = 0.8 + e^{-0.075w(y-1)} + y \qquad (4)$$

$$y = \frac{1}{1 + 10^{-11w^{12}}}, \qquad (5)$$

where $w$ is again the vessel diameter measured in microns, and $H_D$ is the discharge haematocrit. The effective viscosity of the fluid is then given by $\eta = \eta_{\text{vitro}}\eta_0$, where $\eta_0$ is the viscosity of plasma (whose value does not modify the nature of flow solutions).
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Experimental velocity ratio $V_2/V_1$ normalized by that in the reference state, versus inlet haematocrit $H_D$, for different channel widths $w$: (i) thin channels, $w \le 20 \mu m$ and (ii) thicker channels, $w \ge 25 \mu m$. The velocities are measured a few $t_0$ after the flow is established. (b) Theoretical bifurcation diagram for (i) a geometrically symmetric (light) and perturbed network (dark) with $w = 11 \mu m$ and (ii) a geometrically symmetric (light) and perturbed network (dark) with $w = 36 \mu m$. Colored marks correspond to the equilibria supported at $H_D = 0.7$, also reported in Fig. 5(a).</figcaption></figure>

The discharge haematocrit in vessel $i$ is determined by either conservation of RBC flow if vessel $i$ is the outflow of a converging node,

$$Q_i H_i = H_j Q_j + H_k Q_k, \qquad (6)$$

where vessels $j$ and $k$ are the feed vessels, or via plasma skimming if vessel $i$ is the outflow of a diverging node,

$$Q_i H_i = Q_j H_j f(\cdot), \qquad (7)$$

where vessel $j$ is the feed vessel. In general, the plasma skimming function $f$ depends on the flow rate in the daughter vessel relative to the flow in the feed vessel $j$, with lower flows recruiting significantly lower proportions of the total red blood cell balance. Here, we use the empirical law of Pries et al. [22] that also accounts for differences in red cell distribution as it relates to vessel diameters and feed haematocrit:

$$f = \frac{1}{Q} \begin{cases} 0, & Q < Q_0 \\ \left[ 1 + \exp \left( A - B \log \left( \frac{Q-Q_0}{1-Q-Q_0} \right) \right) \right]^{-1}, & Q_0 \le Q \le 1 - Q_0 \\ 1, & Q > 1 - Q_0. \end{cases} \qquad (8)$$

where $Q = |Q_i/Q_j|$ is the fractional flow entering daughter branch $i$ instead of the alternative daughter branch $i'$,
<figure><img src="image_5.png" /><figcaption>FIG. 4. Box plots of number of equilibria in geometrically perturbed networks of reference diameter 11 µm, with diamonds indicating outlier observations. We generated 100 random networks with vessel diameters perturbed by ε up to 10%. For each, we sample 200 different initial conditions and iterate until convergence. An equilibrium is saved only if the solution quality is high and the equilibria are well differentiated from one another.</figcaption></figure>

and the empirical constants have been constructed to capture experimental data:

$$A = \frac{6.96}{w_i} \ln \left( \frac{w_i}{w_{i'}} \right) \qquad (9)$$

$$B = 1 + 6.98 \left( \frac{1 - H_j}{w_j} \right) \qquad (10)$$

$$Q_0 = \frac{0.4}{w_j}. \qquad (11)$$

In our situation where all channels have equal widths, A = 0. Each resistance featured in our three Kirchoff equations can be expressed solely as a function of the three free flows, and so we have closed the system of nonlinear equations. The determination of the equilibrium flows and haematocrits in networks has been well-detailed elsewhere [27], and we use numerical continuation methods to track the equilibrium solutions as we change system parameters. The fully symmetric network exhibits non-generic bifurcations [49], and so we also consider geometrically perturbed systems that more closely match experiments; here we introduce random deviations of less than 1% to all vessel diameters.

We show in Fig. 3(b) the equilibrium solutions for two 2-bridge networks of channel diameters w = 11 and 36 µm respectively, for different inlet haematocrit values. Notice that for the geometrically perturbed system the solution curves become asymmetric but otherwise retain the same features as the fully symmetric case. We have confirmed the generality of these predictions by conducting a statistical study with perturbed nominal resistances via the diameters: $\bar{w} = 11 + u\epsilon$, where u is uniformly sampled from the unit interval and $\epsilon$ is the strength of the perturbation. For each value of $\epsilon$, we generate 100 such nominal resistance profiles. For each nominal resistance profile, we generate a collection of equilibria by generating 200 haematocrit profiles (sampled uniformly randomly from the unit interval for each vessel) and use each of these distributions as the initial condition for a standard nonlinear root finder for the equilibrium equations $F(\mathbf{h}) = \mathbf{0}$, where $\mathbf{F}$ is the equilibrium relation and $\mathbf{h}$ is the vector of haematocrits in each vessel. We keep a candidate $\mathbf{h}^*$ only if the solution quality is high ($\|F(\mathbf{h}^*)\|_2 \le 1 \times 10^{-6}$) and the candidate is well differentiated from all others ($\|\mathbf{h}^* - \mathbf{h}^{*'}\|_2 \ge 1 \times 10^{-3}$). In Fig. 4 we show the results of this analysis. We see that multiple equilibria are supported across all perturbation strengths, but that the number of total equilibria that are supported decreases as the perturbation strength is increased.

The bifurcation diagram of Fig. 3(b) confirms the presence of a threshold in both channel diameter and inlet haematocrit for the existence of non-trivial equilibrium solutions in the network. Note that the threshold haematocrit at which multiple solutions appear in the simulation results are noticeably higher than in experiments ($H_D \simeq 0.62$)
<figure><img src="image_6.png" /><figcaption>FIG. 5. (a): Stability of the equilibria supported by the symmetric networks (semitransparent marks) and perturbed networks (opaque marks) for $H_D = 0.7$ (circular marks) and $H_D = 0.8$ (square marks) with $w = 11 \mu m$ in all cases. The traces represent the trajectories from the time domain simulations in the geometrically perturbed network shown in (b): velocity ratios are plotted as a function of normalized time $\tau = t/t_0$. After ramping up, the inlet haematocrit is $H_D = 0.7$ for $20 < \tau < 50$, then $H_D = 0.8$ for 29 transits (dashed trace) or 33 transits (dot-dashed trace), before returning to $H_D = 0.7$.</figcaption></figure>

vs. $H_D \simeq 0.45$ for $w = 11 \mu m$). This may be related to the choice of the viscosity [50] and plasma skimming models established from rather scattered data obtained in cylindrical channels, that should not be expected to perfectly describe the behavior in channels with rectangular cross-sections. At $H_D = 0.7$ and $w = 11 \mu m$, there are five distinct equilibria supported by the perturbed network, and each of these may be characterized by the flow direction in the bridges — these are shown in the $(V_2/V_1, V_3/V_2)$ space in Fig. 5(a). We can characterize the system by the flow direction in the bridges, e.g. Up-Up (UU), Up-Down (UD), etc.

Stability of the equilibrium flow distribution is determined by linearizing the dynamic advection model for RBC transport through the network,

$$ \frac{\partial H_i}{\partial t} + V_i \frac{\partial H_i}{\partial x_i} = 0, \qquad (12) $$

where $V_i$ is the advection speed through vessel i. The result is a transcendental characteristic equation,

$$ 1 = \sum_j \frac{c_j}{\tau_j} (1 - e^{-\lambda \tau_j}), \qquad (13) $$

where $\tau_j$ is the transit time of every partial flow pathway through the network and $c_j$ is determined by linearizing the effects of plasma skimming and viscosity on the network response [26]. The eigenvalues $\lambda = \sigma + i\omega$ are complex and are challenging to compute in general. For illustration purposes, we show an example in Fig. 6 of the real and imaginary contours of Eq. 13 evaluated at $(V_2/V_1, V_3/V_2) \approx (1.068, 1.090)$ in the perturbed system, i.e., a green mark in Fig. 5. Eigenvalues correspond to the intersection points of the two curve sets, and so we conclude that this equilibrium is a stable spiral, as its eigenvalues with largest real part have $\sigma < 0$ and $\omega \neq 0$. Similarly, an unstable spiral would have dominant eigenvalues with $\sigma > 0$ and $\omega \neq 0$, whereas a saddle would have dominant eigenvalue with $\sigma > 0$ and $\omega = 0$.

We represent in Fig. 5(a) these different types of equilibria in the $(V_2/V_1, V_3/V_2)$ space, for two selected haematocrits and $w = 11 \mu m$. The stable states correspond to UU and DD configurations, which is in agreement with the experimental observation depicted in Fig. 2(b): for the narrowest channel, if we exclude the states where the flux in the bridges are within the experimental uncertainty established while determining the reference state, the probability to be in UU or DD state is 78%, to be compared with a 22% likelihood to be in a UD or DU configuration. The system indeed seems to spend more time in the flow configurations that correspond to the stable states predicted by stability analysis. Generally speaking, the long transit times through the network (due to the slow flow in the bridges) partly dictates the linear response to perturbations and leads to growth (decay) rates that are very large
<figure><img src="image_7.png" /><figcaption>FIG. 6. Zero contours of the real (red traces) and imaginary (blue traces) parts of Eq. 13 evaluated at $(V_2/V_1, V_3/V_2) \approx (1.072, 1.083)$ in the perturbed system of Fig. 3(b-i). Black marks indicate eigenvalues of the system.</figcaption></figure>

compared to the bulk transit time in the network. We therefore expect to find slowly evolving dynamics coupled with fast transitions mediated by the saddle points.

We confirm these theoretical predictions with independent time-domain simulations of the full nonlinear advection equations, as shown for the perturbed network in Fig. 5(b). The inlet haematocrit is initially ramped-up from $H_D = 0$ to $H_D = 0.7$ over the course of $\tau = t/t_0 = 20$ transit times and the system converges to the stable DD-state. At $\tau = 50$ we instantaneously increase the inlet haematocrit to $H_D = 0.8$ and the system gets kicked toward a sustained oscillation in the vicinity of an unstable spiral. If we return to $H_D = 0.7$ after 29 transit times, the network is attracted to the stable UU-state. In contrast, if we return to $H_D = 0.7$ after 33 transit times, the network is attracted to the stable DD-state, i.e. the outcome depends on the interplay between the basin of attraction of each stable equilibrium and the state of the system when the inlet haematocrit is changed.

The time domain simulation shown in Fig. 5(b) exhibits transitions between states that resemble that seen in experiments (Fig. 2(a)) and suggests that the system may converge, after many bulk transit times, to a stable, steady state or a sustained oscillation.

## RELEVANCE TO PHYSIOLOGICAL FLOWS

Importantly, the values of haematocrit and channel sizes at which symmetry-breaking occurs in our experiments are near the range of human microcirculatory physiological conditions ($H_D \approx 0.5$ and $w \approx 10$ µm, corresponding to an average tube haematocrit $H_T \approx 0.37$). The discharge haematocrit $H_D$ is a conserved flux ratio, indicating that microvessels are typically perfused with this haematocrit, with large fluctuations around this average value due to phase separation [51]. This makes this study relevant for a better understanding of RBC distribution heterogeneities in *in vivo* networks.

Note that while we initially used in our simulations an *in vitro* viscosity law [50] to tentatively match *in vitro* experiments, similar computations using *in vivo* data laws show a lower threshold within physiological ranges: Pries et al. have also supplied a widely used quantification of the Fåhræus-Lindqvist effect *in vivo* [50]:

$$ \eta_{\text{vivo}} = \left[ 1 + z(\eta^* - 1) \frac{(1 - H^D)^C - 1}{(1 - 0.45)^C - 1} \right] z \quad (14) $$

$$ \eta^* = 6e^{-0.085w} + 3.2 - 2.44e^{-0.06w^{0.645}} \quad (15) $$

$$ z = \left( \frac{w}{w - 1.1} \right)^2, \quad (16) $$

where C and y are identical to the *in vitro* formulation, and w and $H_D$ are again the vessel diameter in microns
<figure><img src="image_8.png" /><figcaption>FIG. 7. Bifurcation diagram for the symmetric network (semi-transparent marks) and geometrically perturbed network (opaque marks) with $w = 11$ µm and *in vivo* viscosity.</figcaption></figure>

and the discharge haematocrit, respectively. Due to this change in viscosity law, and more precisely to a larger logarithmic derivative [52, 53], we note a significantly lower multistability threshold, $H_D \approx 0.41$.

The existence of a threshold well within the physiological range is consistent with previous studies on complex networks [28, 54]. It also suggests the existence of lower thresholds *in vivo* compared to that we exhibited in our *in vitro* experiments.

DISCUSSION AND CONCLUSION

By experimentally characterizing the evolution of velocity and haematocrit distributions in a simple, nearly symmetric network made of two bridges, we demonstrate that multiple, non-trivial flow solutions emerge and that haematocrit and channel sizes are the main parameters triggering this emergence, with thresholds in a physiologically relevant range. These multiple solutions are favored by the existence of multiple paths: blood flow in a single loop with no bridge exhibits a single solution [47]. Of critical importance, these solutions are predicted by a continuous model based on few ingredients without needing to invoke discrete events as was the case in several previous studies [5, 6, 32]. This suggests to use perturbation methods of this continuous model to account for the flow fluctuations in microcirculation.

While in principle the system should converge to a steady state or a sustained oscillation, we have also highlighted that the long-term outcome, and the values of the associated flux, are highly sensitive to both the haematocrit and geometrical imperfections, as exemplified by the difference between the symmetric and geometrically perturbed networks in Fig. 5(a), for which we simply introduced a diameter variation of 1%. These imperfections, or haematocrit fluctuations, can also be dynamically generated by local clogging/adhesion of circulating cells and should be considered as always present, yet fluctuating, in a network. Since at the same time the system always converges slowly due to the control of the dynamics by lateral bridges and their low velocity flux, we therefore expect that such a system will mostly be in a non-steady situation. These findings suggest that the usual assumption of a single stable state would deserve to be revisited, even in rather simple networks, and that functional "disordered states" may emerge even in the absence of pathologies.

In our simplified mesh-like network, the identified solutions favor RBC flow in the branches that do not belong to the most direct path and help homogenize the distribution of cells, especially when it fluctuates between different configurations. In the meantime, the cells flowing in these transverse paths will spend more time in the network. We therefore expect the distribution of transit times to broaden towards large values. The potential analogy with the broadening of distribution of transit times obtained as multiple paths in space are explored in mesh-like networks [12] deserves to be explored.

MA, GC and TP acknowledge financial support by LabEx Tec 21 (Investissements d'Avenir - grant agreement ANR-11-LABX-0030). GC and TP are grateful to B. de Vicente, G. Lento and M. Leonardo for preliminary experiments,
and to M. van Melle-Gateau for technical support for photolithography. TP and GC acknowledge support from CNES
(Centre National d'Etudes Spatiales, DAR ID 4759).

Experimental Methods

Sample preparation

Blood samples were collected at the *Etablissement Français du Sang* (French Blood Agency, EFS) and stored in a citrate solution to prevent coagulation. The blood preparation protocol follows the prescriptions of Merlo et al. [41]. Cells are washed 3 times in Phosphate-Buffered Saline (PBS, P4417 from Sigma Aldrich) by centrifugation, before being mixed with a buffer solution made of 68.5% PBS, 31.5% Optiprep (Axis Shield) and 2 g/L Bovine Serum Albumin (BSA, A7906 from Sigma). The haematocrit (volume fraction in RBCs) $H_{D0}$ of the preparation to be injected in the network is computed from the initial haematocrit $H_i$ in collected tubes provided by EFS and the initial sample volume $V_i$ and final prepared volume $V_0$ through: $H_{D0} = H_iV_i/V_0$. Samples are used straight after preparation.

Experimental set-up

The chip is fabricated by pouring and curing polydimethylsiloxane (PDMS) on a mould produced by direct soft
lithography (Kloé, Dilase 250 and photoresist SU8 1070, Gersteltec). It is then bonded on a glass slide after plasma
treatment. All channels in the 2-bridge network have equal length L = 0.5 mm, equal height h = 30 µm, and equal
width w that was varied from 11 to 36 µm.

Fig. 1(a) shows the 2-bridge network and upstream channels with two distinct inlets fed respectively either by highly concentrated RBCs (haematocrit $H_{D0} \approx 60-70\%$), or buffer solution. They merge in a converging bifurcation located 11 mm upstream the network ($\approx 300$ times the width of the channel), ensuring an homogeneous mixing of the RBCs with the buffer solution in the $\approx 38 \times 30 \mu m$ cross section serpentine, thanks to both wall-induced lift and shear-induced diffusion [16]. The channels between the inlets and this bifurcation have a length of 1.4 mm, width $38 \mu m$ and thickness $h = 30 \mu m$. The flow rate in each inlet is varied using a pressure controller *Elveflow OB1* offering a good stability over time in the pressure range 0 – 2 bar.

The first step consists in filling the network with the buffer solution for 2 hours, allowing for adsorption of BSA on the walls of the channels to prevent cell adhesion. Then both inlet are pressurized at the desired values in order to reach the desired velocities and the targeted haematocrit $H_D$ to feed the 2-bridge network. These values are typically, for the RBCs inlet and buffer inlet respectively: 30 mbar - 30 mbar for intermediate haematocrit (as in Fig. 1(a)), 45 mbar - 15 mbar to reach $Q_{buffer} = 0$ such that $H_D = H_{D0}$ ((as in Fig. 1(b)), 15 mbar - 45 mbar to reach $Q_{RBCs} = 0$ for rinsing procedures. The sum of the inlet pressures is kept constant so as to ensure that the total flow rate stays in a narrow range, therefore the RBC maximal velocities in the network are always in the physiological range 0.5-2 mm.

Fundamental aspects of this set-up are first that different and continuously adjustable haematocrits can be injected
in the network without preparing and handling different samples in the reservoirs, thus allowing for an easy control
of the fundamental parameter $H_D$. Second, rinsing of the network is achieved in a smooth way and in a reasonable
time of a few seconds without increasing hydrodynamic stresses in the network. The ability to easily rinse the system
is important since observed flow patterns could be the consequence of undesired adhesion events of cells, or clogging
by dusts, which cannot be detected when the network is filled with a highly concentrated suspension. Therefore,
after each measurement, rinsing is systematically made by setting $Q_{RBCs} = 0$ while keeping the sum of pressures at
the inlets constant in order not to create additional stress that would remove the clogging without noticing it. This
enables a visual check of possible cell adhesion or clogging in the network. In that case, the whole data set acquired
from the last rinsing to this point is considered as untrustful and discarded.

For both short and long time measurements, we set the desired ratio of RBCs to buffer flow rates after rinsing, then
wait around 30 s before triggering the acquisition. This time is generally enough to reach a quite constant haematocrit
at the inlet of the network inlet (see Fig. 9).
## Image acquisition

The microfluidic chip is set on the stage of an Olympus IX71 inverted microscope, illuminated with an external light source. The image is magnified 32 times, resulting in a pixel size of 1.08 µm. For each acquisition, i.e. each measurement of velocities in the network, 260 images are recorded using a Phantom V2511 camera of 1280×800 pixels resolution, an exposure time $t_{exp} = 5$ µm and a sampling frequency $f_{ech} = 1000$ Hz. We ensured that the number of images acquired at this frequency is enough to get converged data for the velocity. For long measurements, an Arduino microcontroller triggers the camera every 5 seconds. Due to RAM limitations, this time step is increased to 20 seconds after ≈ 400 seconds of acquisition.

### Image processing: haematocrit

Two types of haematocrit are usually defined. The tube haematocrit corresponds to the instantaneous volume fraction of RBCs averaged over the channel cross-section A:

$$H_T = \frac{1}{A} \int HdA, \qquad (17)$$

where H is the local volume fraction. The discharge haematocrit, also called reservoir haematocrit, is the volume fraction of RBCs weighted by their velocity:

$$H_D = \frac{1}{Q} \int uHdA, \qquad (18)$$

where u is the local velocity and $Q = \int u dA$ the global flow rate. $H_D$ is therefore the relative RBC flux.

In a large channel, where the surrounding fluid and the RBCs flow on average at the same velocity, these two parameters coincide. In narrow channels, RBCs tend to accumulate in the center and therefore flow faster on average than the suspending medium. As a consequence, $H_T < H_D$. This effect is called the Fåhræus effect. Contrary to the tube haematocrit, the discharge haematocrit, which is the ratio of the flux of RBCs against the total flux, is a conserved quantity along channel, whatever its section. It is therefore a relevant quantity to be used as a control parameter for comparing 2-bridge networks of different widths w.

It is computed in the inlet area upstream the network (see Fig. 1(a)) by measuring light absorption following a Beer-Lambert approach:

The mean tube haematocrit is first computed from grey level intensities measured in the region of interest (ROI) whose top and bottom boundaries are set as close as possible from the channel walls. The mean intensity I is obtained by averaging in space and time, in the ROI. The optical density is then computed as:

$$OD = -\log \left[ \frac{I}{\alpha I_r} \right], \qquad (19)$$

where $I_r$ is the reference averaged intensity in the same ROI when no RBCs flow in, obtained from a different acquisition. Since several minutes elapsed between the acquisitions from which I and $I_r$ are computed, slight ambient light variations may occur and should be taken into account. Hence the coefficient $\alpha = I_{PDMS}/I_r$, PDMS in Eq. 19, which is the ratio of the averaged intensity in the PDMS (far away from the channel) computed from the considered acquisition, to that in the reference acquisition at the same location. Eventually, assuming a Beer-Lambert law as in [16, 43, 55, 56], the tube haematocrit in the ROI is computed as follows:

$$H_T = H_{T0} \times OD/OD_0, \qquad (20)$$

where OD₀ = optical density obtained when non diluted suspension is injected ($Q_{buffer} = 0$): the tube haematocrit is then $H_{T0}$, corresponding to the known discharge haematocrit $H_{D0}$.

The conversion between tube and discharge haematocrit can then be achieved using an empirical relation established by Pries [19]:

$$H_T = H_D(H_D + (1 - H_D)X), \qquad (21)$$
<figure><img src="image_9.png" /><figcaption>FIG. 8. Illustration of the method for fluid velocimetry. $\Delta z$ is the distance between the two windows of width $w_s$, and the time step between the two images is of the order of 1 ms.</figcaption></figure>

where $X = 1 + 1.7e^{-0.35d} - 0.6e^{-0.01d}$, and $d$ is the diameter of a cylindrical channel, measured in microns, for which this law has been proposed initially. It has been shown in [57] that viscosity laws in rectangular cross section channels, that are also a consequence of cell centering, can be well described by the law also proposed by Pries et al. [19] for cylindrical channels, by setting $d = \min(w, h)$. Following this idea, we use Eq. 21 with $d = h$.

Finally, $H_{T0}$ is obtained from $H_{D0}$ through Eq. 21, and $H_D$ is obtained from $H_T$, which was calculated thanks to Eq. 20, by solving Eq. 21, which provides

$$H_D = -\frac{X}{2-2X} + \left[ \left( \frac{X}{2-2X} \right)^2 + \frac{H_T}{1-X} \right]^{1/2}. \qquad (22)$$

### Image processing: fluid velocity

The velocimetry consists in calculating an averaged point to point space-correlation of intensities measured at two locations spaced by a varying distance $\Delta z$ along the direction of the flow, separated by a time step $\Delta t$ (see Fig. 8). By integrating the space-correlation function over windows of length $w_s$ and width $w$, and considering all the pairs of images available in the recorded sequence, one eventually gets a correlation function whose maximum corresponds to the maximal fluid velocity within the interrogation window.

For a given pair of images $n, n+m$ corresponding to times $t = n/f_{ech}$, $t + \Delta t = (n+m)/f_{ech}$ ($f_{ech}$ is the sampling frequency), the measured signals are the grey levels of the pixels located at $(z, y)$ and $(z + \Delta z, y)$ respectively:

$$\begin{aligned} s_1 &= I(z, y, n), \\ s_2 &= I(z + \Delta z, y, n+m) \end{aligned} \qquad (23)$$

for all the pixels of the interrogation window on each image. The averaged normalized correlation between the two signals is computed on the full range of images and over the areas of analysis:

$$C(\Delta z) = \frac{\sum_{n,y,z} s_1 s_2}{\sqrt{\sum_{n,y,z} s_1^2} \sqrt{\sum_{n,y,z} s_2^2}}, \qquad (24)$$

where $N$ is the number of images of the acquisition, $\sum_{n,y,z} \equiv \sum_{n=1}^{N-m} \sum_{y=y_0}^{y_0+w} \sum_{z=z_0}^{z_0+w_s}$ denotes the sum over the available pairs of images and over the whole interrogation window, and $(y_0, z_0)$ is the location of the bottom-left corner of the interrogation window on the first image of each pair.

$C_{\Delta z}$ is computed for several values of $\Delta z$, varied from -20 to +20 pixels, and the fluid velocity $V_f$ is deduced from the value of $\Delta z_0$ that corresponds to the position of the maximum of the correlation function (obtained by a local Gaussian fit to reach subpixel precision):

$$V_f = \frac{\Delta z_0 \delta}{\Delta t}, \qquad (25)$$
<figure><img src="image_10.png" /><figcaption>FIG. 9. (a) Inlet haematocrits and velocities as a function of time, for the data shown in Fig. 2(a). (b) Experimental velocity ratios $V_2/V_1$ and $V_3/V_2$ as a function of normalized time, for a channel width $w$ of 11 µm. Both graphs correspond to two additional experimental realizations with similar parameters as in Fig. 2(a). The inlet haematocrit $H_D$ is constant in time ($H_D \simeq 0.65 \pm 1\%$). Relative concentrations in each bridge are shown in gray scale normalized by the maximum concentration measured in the bridge. Colored bands identify the reference state, accounting for measurement fluctuations.</figcaption></figure>

where $\delta$ is the pixel size. As in [58] where a time-correlation method was used, we checked that the measured velocity corresponds to the maximal velocity in the $xy$ plane (i.e. central velocity) by applying our method to simulated Poiseuille flows of objects of elliptic cross section with random position and orientations.

In practice, the time interval between the two images must be small enough to get a strong correlation between the signals measured in the two windows, but large enough so that the relative error on the detection of the position of the maximum is small. In practice, this would be reached when the cells are displaced by a distance of the order of a RBC diameter between the two acquisition times, after what correlation is lost because of the complexity of cell dynamics in concentrated suspensions. With $\Delta z_0$ of the order 10 µm, a velocity of the order 1 µm, $\Delta t$ should be around 10 ms. Finally, a quick sensitivity analysis has been performed in the range $2 \le \Delta t \le 5$ ms for a given dataset: the computed velocity varies within $\pm 5\%$. Therefore, $\Delta t$ is initially set at 3 ms. If the top of the correlation function cannot be fitted by a Gaussian function (i.e. the correlation function does not show a clear peak), the process is repeated by increasing $\Delta t$ by steps of 1 ms, and up to 10 ms. Furthermore, the process is also sensitive to the number of pairs of images processed. We thus studied the convergence of the computed velocity for $N-m=30, 60$ and $120$ pairs, and observed that $N-m \ge 60$ is enough to get converged data.

## Reference state determination

In order to compare networks of different widths to build up Fig. 3(a), we determine a "reference state" that would be the symmetric state (no flow in the bridges) for a perfect network and is in practice slightly asymmetric because of the imperfections due to the microfabrication process (precision of photolithography, small deformations when bonding the PDMS chip to the glass slide,...).

When $H_D$ was varied, as for the state diagram (short-time measurements) of Fig. 3(a), we simply define the
reference velocity ratio as the mean values when $H_D \le 0.3$, i.e when the fluid is expected to behave as a simple, Newtonian, fluid. However, when processing time series (Fig. 2), in the absence of low haematocrit data, we first looked for data filling the condition $H_{bridge}/\max(H_{bridge}) < 0.05H_{D0}$, i.e. situations with almost no cells in the bridges. If no data could satisfy this condition, which can happen if the bridges are initially filled with cells that then do not move, we then searched data verifying $[H_{bridge}/\max(H_{bridge}) < 0.6H_{D0} \& 0.9 \le V_{n+1}/V_n \le 1.1]$, i.e. situations with small enough flow in the bridges. Taking the mean values for the selected data, we obtain the reference values for the velocity ratios. The colored intervals in Fig. 2 in the main paper correspond to these values $\pm$ the standard deviation. For experiments where low $H_D$ data were available, we found that this method induces a deviation of no more than 8% compared to the expected value.

Additional data for long time measurements

We plot in Fig. 9 the inlet velocities and haematocrits as a function of time for the experiments leading to the data of Fig. 2(a); these values converge to a quasi steady one in about 100 s. The data shown in Fig. 2(a) are taken when this plateau is reached, i.e when variations do not exceed 10 % of the mean value.

Fig. 9 shows two other long-time measurements in the $w = 11~\mu\text{m}$ network. They complement the graph shown in Fig. 2(a) in the main paper.

* thomas.podgorski@univ-grenoble-alpes.fr

† gwennou.coupier@univ-grenoble-alpes.fr

[1] M. Ruiz-Garcia and E. Katifori, “Emergent dynamics in excitable flow systems,” Physical Review E 103, 062301 (2021).

[2] C. F. Daganzo, V. V. Gayah, and E. J. Gonzales, “Macroscopic relations of urban traffic variables: Bifurcations, multi-valuedness and instability,” Transportation Research Part B: Methodological 45, 278–288 (2011).

[3] M. Prakash and N. Gershenfeld, “Microfluidic bubble logic,” Science 315, 832 (2007).

[4] M. Schindler and A. Ajdari, “Droplet traffic in microfluidic networks: A simple model for understanding and designing,” Phys. Rev. Lett. 100, 044501 (2008).

[5] O. Cybulski, S. Jakiela, and P. Garstecki, “Between giant oscillations and uniform distribution of droplets: The role of varying lumen of channels in microfluidic networks,” Phys. Rev. E 92, 063008 (2015).

[6] O. Cybulski, P. Garstecki, and B.A. Grzybowski, “Oscillating droplet trains in microfluidic networks and their suppression in blood flow,” Nat. Phys. 15, 706–713 (2019).

[7] M. F. Kiani, A. R. Pries, L. L. Hsu, I. H. Sarelius, and G. R. Cokelet, “Fluctuations in microvascular blood flow parameters caused by hemodynamic mechanisms,” American Journal of Physiology-Heart and Circulatory Physiology 266, H1822–H1828 (1994).

[8] L. V. Mezentseva, S. S. Pertsov, and V. K. Hugaeva, “A comparative analysis of the persistence of capillary blood flow oscillations in the left and right rat kidneys,” Biophysics 61, 656–660 (2016).

[9] P. Bagher and S. S. Segal, “Regulation of blood flow in the microcirculation: role of conducted vasodilation,” Acta Physiologica 202, 271–284 (2011).

[10] A. S. Popel and P. C. Johnson, “Microcirculation and hemorheology,” Annu. Rev. Fluid Mech. 37, 43–69 (2005).

[11] Sylvie Lorthois and Francis Cassot, “Fractal analysis of vascular networks: Insights from morphogenesis,” Journal of Theoretical Biology 262, 614–633 (2010).

[12] F. Goirand, T. Le Borgne, and S. Lorthois, “Network-driven anomalous transport is a fundamental component of brain microvascular dysfunction,” Nat. Commun. 12, 7295 (2021).

[13] J. Dupire, M. Socol, and A. Viallat, “Full dynamics of a red blood cell in shear flow,” Proc. Nat. Acad. Sci. USA 109, 20808 (2012).

[14] J. Mauer, S. Mendez, L. Lanotte, F. Nicoud, M. Abkarian, G. Gompper, and D. A. Fedosov, “Flow-induced transitions of red blood cell shapes under shear,” Phys. Rev. Lett. 121, 118103 (2018).

[15] C. Minetti, V. Audemar, T. Podgorski, and G. Coupier, “Dynamics of a large population of red blood cells under shear flow,” J. Fluid Mech. 864, 408 (2019).

[16] X. Grandchamp, G. Coupier, A. Srivastav, C. Minetti, and T. Podgorski, “Lift and down-gradient shear-induced diffusion in red blood cell suspensions,” Phys. Rev. Lett. 110, 108101 (2013).

[17] M. Brust, O. Aouane, M. Thiébaud, D. Flormann, C. Verdier, L. Kaestner, M. Laschke, H. Selmi, A. Benyoussef, T. Podgorski, G. Coupier, C. Misbah, and C. Wagner, “The plasma protein fibrinogen stabilizes clusters of red blood cells in microcapillary flows,” Sci. Rep. 4, 4348 (2014).

[18] Walter H. Reinhart, Nathaniel Z. Piety, and Sergey S. Shevkoplyas, “Influence of red blood cell aggregation on perfusion of an artificial microvascular network,” Microcirculation 24, e12317 (2017).

[19] A. Pries, N. Neuhaus, and P. Gaehtgens, “Blood viscosity in tube flow: dependence on diameter and hematocrit,” Am. J. Physiol. 20, H1770–H1778 (1992).
[20] J. W. Dellimore, M. J. Dunlop, and P. B. Canham, "Ratio of cells and plasma in blood flowing past branches in small plastic channels," Am. J. Physiol. Heart Circ. Physiol. **244**, H635-H643 (1983).

[21] B. M. Fenton, R. T. Carr, and G. R. Cokelet, "Nonuniform red cell distribution in 20 to 100 µm bifurcations," Microvasc. Res. **29**, 103-126 (1985).

[22] A. R. Pries, K. Ley, M. Claassen, and P. Gaethgens, "Red cell distribution at microvascular bifurcations," Microvasc. Res. **38**, 81-101 (1989).

[23] V. Doyeux, T. Podgorski, S. Peponas, M. Ismail, and G. Coupier, "Spheres in the vicinity of a bifurcation: elucidating the Zweifach-Fung effect," J. Fluid Mech. **674**, 359 (2011).

[24] P. Balogh and P. Bagchi, "Analysis of red blood cell partitioning at bifurcations in simulated microvascular networks," Phys. Fluids **30**, 051902 (2018).

[25] J. Geddes, R. Carr, N.J. Karst, and F. Wu, "The onset of oscillations in microvascular blood flow," SIAM J. Appl. Dyn. Syst. **6**, 694-727 (2007).

[26] J.M. Davis and C. Pozrikidis, "Self-sustained oscillations in blood flow through a honeycomb capillary network," Bull. Math. Biol. **76**, 2217-2237 (2014).

[27] N.J. Karst, B.D. Storey, and J.B. Geddes, "Oscillations and multiple equilibria in microvascular blood flow," Bull. Math. Biol. **77**, 1377-1400 (2015).

[28] Nathaniel J. Karst, John B. Geddes, and Russell T. Carr, "Model microvascular networks can have many equilibria," Bull. Math. Biol. **79**, 662-681 (2017).

[29] J. M. Davis, "On the linear stability of blood flow through model capillary networks," Bull. Math. Biol. **76**, 2985-3015 (2014).

[30] Yaron Ben-Ami, George W Atkinson, Joe M Pitt-Francis, Philip K Maini, and Helen M Byrne, "Structural features of microvascular networks trigger blood flow oscillations," Bulletin of Mathematical Biology **84**, 85 (2022).

[31] Noemie Boissier, Dirk Drasdo, and Irene E. Vignon-Clementel, "Simulation of a detoxifying organ function: Focus on hemodynamics modeling and convection-reaction numerical simulation in microcirculatory networks," Int. J. Num. Meth. Biomed. Eng. **37**, e3422 (2021).

[32] Guansheng Li, Ting Ye, Zehong Xia, Sitong Wang, and Ziwei Zhu, "Analysis and prediction of hematocrit in microvascular networks," International Journal of Engineering Science **191**, 103901 (2023).

[33] Angela d'Esposito, Paul W. Sweeney, Morium Ali, Magdy Saleh, Rajiv Ramasawmy, Thomas A. Roberts, Giulia Agliardi, Adrien Desjardins, Mark F. Lythgoe, R. Barbara Pedley, Rebecca Shipley, and Simon Walker-Samuel, "Computational fluid dynamics with imaging of cleared tissue and of in vivo perfusion predicts drug uptake and treatment responses in tumours," Nature Biomedical Engineering **2**, 773-787 (2018).

[34] Franca Schmid, Matthew J. P. Barrett, Dominik Obrist, Bruno Weber, and Patrick Jenny, "Red blood cells stabilize flow in brain microvascular networks," PLOS Computational Biology **15**, 1-29 (2019).

[35] A. Mantegazza, F. Clavica, and D. Obrist, "In vitro investigations of red blood cell phase separation in a complex microchannel network," Biomicrofluidics **14**, 014101 (2020).

[36] Miguel O. Bernabeu, Jakub Köry, James A. Grogan, Bostjan Markelc, Albert Beardo, Mayeul d'Avezac, Romain Enjalbert, Jakob Kaeppler, Nicholas Daly, James Hetherington, Timm Krüger, Philip K. Maini, Joe M. Pitt-Francis, Ruth J. Muschel, Tomás Alarcón, and Helen M. Byrne, "Abnormal morphology biases hematocrit distribution in tumor vasculature and contributes to heterogeneity in tissue oxygenation," Proc. Nat. Acad. Sci. **117**, 27811-27819 (2020).

[37] R. Enjalbert, T. Krüger, and M. O. Bernabeu, "Effect of vessel compression on blood flow in microvascular networks and its implications for tumour tissue hypoxia," Commun. Phys. **7**, 49 (2024).

[38] O. Forouzan, X. Yang, J. M. Sosa, J. M. Burns, and S. S. Shevkoplyas, "Spontaneous oscillations of capillary blood flow in artificial microvascular networks," Microvasc. Res. **84**, 123-132 (2012).

[39] Shyr-Shea Chang, Shenyinying Tu, Kyung In Baek, Andrew Pietersen, Yu-Hsiu Liu, Van M. Savage, Sheng-Ping L. Hwang, Tzung K. Hsiai, and Marcus Roper, "Optimal occlusion uniformly partitions red blood cells fluxes within a microvascular network," PLOS Computational Biology **13**, 1-22 (2017).

[40] Sergey S. Shevkoplyas, Tatsuro Yoshida, Sean C. Gifford, and Mark W. Bitensky, "Direct measurement of the impact of impaired erythrocyte deformability on microvascular network perfusion in a microfluidic device," Lab Chip **6**, 914-920 (2006).

[41] Adlan Merlo, Sylvain Losserand, F. Yaya, Philippe Connes, Magalie Faivre, Sylvie Lorthois, Christophe Minetti, Elie Nader, Thomas Podgorski, Céline Renoux, Gwennou Coupier, and Emilie Franceschini, "Influence of storage and buffer composition on the mechanical behavior of flowing red blood cells," Biophysical Journal **122**, 360-373 (2023).

[42] See supplemental material at [to be inserted by editor].

[43] A. Merlo, M. Berg, P. Duru, F. Risso, Y. Davit, and S. Lorthois, "A few upstream bifurcations drive the spatial distribution of red blood cells in model microfluidic networks," Soft Matter **18**, 1463-1478 (2022).

[44] R. Fahraeus, "The suspension stability of the blood," Physiological Reviews **IX**, 241-274 (1929).

[45] D. A. Fedosov, B. Caswell, A. S. Popel, and G. E. Karniadakis, "Blood flow and cell-free layer in microvessels," Microcirculation **17**, 615-628 (2010).

[46] S. Losserand, G. Coupier, and T. Podgorski, "Migration velocity of red blood cells in microchannels," Microvasc. Res. **124**, 30 (2019).

[47] Oliver Jensen, Serban Pop, Sarah Waters, and Giles Richardson, "Shock formation and nonlinear dispersion in a microvascular network," Mathematical Medicine and Biology **24**, 379-400 (2007).

[48] Y. Ben-Ami, G. W. Atkinson, J. M. Pitt-Francis, P. K. Maini, and H. M. Byrne, "Structural features of microvascular networks trigger blood flow oscillations," Bulletin of Mathematical Biology **84**, 85 (2022).
[49] Y Kuznetsov, Elements of Applied Bifurcation Theory (Springer Verlag, 1998).

[50] A. R. Pries, T. W. Secomb, T. Gessner, M. B. Sperandio, J. F. Gross, and P. Gaehtgens, "Resistance to blood flow in microvessels in vivo." Circ. Res. **75**, 904–915 (1994).

[51] O.K. Baskurt, M.R. Hardeman, M.W. Rampling, and H.J. Meiselman, *Handbook of hemorheology and hemodynamics* (IOS Press, 2007).

[52] J. B. Geddes, B. D. Storey, D. Gardner, and R. T. Carr, "Bistability in a simple fluid network due to viscosity contrast," Phys. Rev. E **81**, 046316 (2010).

[53] Casey M. Karst, Brian D. Storey, and John B. Geddes, "Laminar flow of two miscible fluids in a simple network," Phys. Fluids **25**, 033601 (2013).

[54] Nathaniel J. Karst and John B. Geddes, "Modeling transit time distributions in microvascular networks," J. Theor. Biol. **572**, 111584 (2023).

[55] J. M. Sherwood, E. Kaliviotis, J. Dusting, and S. Balabani, "Hematocrit, viscosity and velocity distributions of aggregating and non-aggregating blood in bifurcating microchannel," *Biomicrofluidics* **6**, 024119 (2012).

[56] S. Roman, A. Merlo, P. Duru, F. Risso, and S. Lorthois, "Going beyond 20 µm-sized channels for studying red blood cell phase separation in microfluidic bifurcations," *Biomicrofluidics* **10**, 034103 (2016).

[57] V. Audemar, T. Podgorski, and G. Coupier, "Rheology and structure of a suspension of deformable particles in plane poiseuille flow," Phys. Fluids **34**, 042013 (2022).

[58] S. Roman, S. Lorthois, P. Duru, and F. Risso, "Velocimetry of red blood cells in microvessels by the dual-slit method: Effect of velocity gradients," *Microvasc. Res.* **84**, 249 (2012).