Structural Features of Microvascular Networks Trigger Blood Flow Oscillations

We analyse mathematical models in order to understand how microstructural features of vascular networks may affect blood flow dynamics, and to identify particular characteristics that promote the onset of self-sustained oscillations. By focusing on a simple three-node motif, we predict that network “redundancy”, in the form of a redundant vessel connecting two main flow-branches, together with differences in haemodynamic resistance in the branches, can promote the emergence of oscillatory dynamics. We use existing mathematical descriptions for blood rheology and haematocrit splitting at vessel branch-points to construct our flow model; we combine numerical simulations and stability analysis to study the dynamics of the three-node network and its relation to the system’s multiple steady-state solutions. While, for the case of equal inlet-pressure conditions, a “trivial” equilibrium solution with no flow in the redundant vessel always exists, we find that it is not stable when other, stable, steady-state attractors exist. In turn, these “nontrivial” steady-state solutions may undergo a Hopf bifurcation into an oscillatory state. We use the branch diameter ratio, together with the inlet haematocrit rate, to construct a two-parameter stability diagram that delineates regimes in which such oscillatory dynamics exist. We show that flow oscillations in this network geometry are only possible when the branch diameters are sufficiently different to allow for a sufficiently large flow in the redundant vessel, which acts as the driving force of the oscillations. These microstructural properties, which were found to promote oscillatory dynamics, could be used to explore sources of flow instability in biological microvascular networks.


Introduction
Sustained oscillations in the microcirculation have been known to occur for some time. They have been observed in vivo (Kimura et al. 1996) and in vitro (Forouzan et al. 2012), and studied via mathematical modelling (for example, Kiani et al. 1994;Carr and Lacoin 2000;Geddes et al. 2010). Such theoretical studies have confirmed the existence of oscillatory solutions that are self-induced, i.e. they emerge in the absence of external forcing or the imposition of oscillatory boundary conditions. In this work, we seek to identify the microstructural features of vascular networks which promote oscillatory dynamics. Understanding such relationships is important in order to predict the behaviour of large-scale vascular networks and the tissue oxygenation they provide. While several theoretical studies have demonstrated that network geometry can affect the emergence of blood flow oscillations (see, for example, Geddes et al. 2007; Davis and Pozrikidis 2014a;Karst et al. 2015), to our knowledge no explicit mechanism has been proposed for how specific structural features of microcapillary networks, together with the inherent properties of blood flow, act to generate oscillatory dynamics. Therefore, in this work, we consider a simple three-node network motif as a model case study to investigate how geometrical features of the network promote oscillatory instability of steady blood flow. Our findings could be used in future work to identify motifs in larger networks that act as sources of instability and trigger oscillatory dynamics.
When blood flows in microcapillaries its viscosity is governed by the concentration of red blood cells (RBCs), such that the hydraulic resistance depends nonlinearly on the haematocrit level, a phenomenon known as the Fåhraeus-Lindqvist effect (Fåhraeus and Lindqvist 1931). On the other hand, the splitting of haematocrit at a vessel branch point depends (nonlinearly) on the partitioning of blood flow between daughter vessels, a phenomenon known as "plasma skimming" (Krough 1921). Together, these two effects result in coupled nonlinear relations between haematocrit concentrations and flow rates in the different vessels of the network. Such intrinsic nonlinearities have been shown to drive the emergence of multiple equilibria and oscillatory dynamics .
The complex rheology of microcapillary blood flow and its strong dependence on the haematocrit and vessel diameter were originally studied by Fåhraeus and Lindqvist (1931). Later on, Pries et al. (1994) derived a widely used mathematical model to quantify these effects. While the Pries et al. (1994) viscosity model is a fundamental element of almost every study of blood flow in the microcirculation, a variety of models for plasma skimming have been used. The haematocrit splitting rules vary from simple, single-parameter equations (Klitzman and Johnson 1982;Fenton et al. 1985) to complicated models based on experimental measurements (Pries et al. 1989) and discrete-RBC simulations (Bernabeu et al. 2020).
The functional form of the haematocrit splitting models has been shown to significantly affect the emergence of self-sustained oscillations; specifically, a dominant factor is the rate at which the haematocrit flux into a specific daughter branch increases as the total flow rate into that branch increases. For example, Pozrikidis (2011, 2014a, b) showed that regular networks, such as honeycomb or tree networks, are only prone to oscillations when physically unrealistic splitting rules, with very large haematocrit-flux gradients, are used. Nevertheless, other research Gardner et al. 2010) has shown that when the more biologically-sound haematocrit splitting model of Pries et al. (1989) is used, multiple equilibria and oscillations may occur if some "redundancy" is introduced into the system; this was achieved by adding a vessel to connect two main flow branches. The term "redundancy" is used here because, in certain conditions, the additional connecting vessel can support a "zero-flow" solution where it transports neither plasma nor haematocrit. These findings suggest that self-induced oscillations are more likely to occur in vascular networks with irregular topological structures, which are, in fact, also characteristic of tumour vasculature (Jain 2005).
Cancer cells influence and respond to local environmental conditions; this leads to rapid and localised angiogenesis, and the formation of networks whose morphologies differ dramatically from those of healthy tissues (Jain 2005). It is hypothesised that the abnormal structure of tumour vasculature leads to spatio-temporal variations in blood flow and haematocrit distribution, which manifest at the macroscopic level as cycling hypoxia (Michiels et al. 2016;Gillies et al. 2018). This phenomenon is characterized by periodic episodes of oxygen deprivation, followed by periods of reoxygenation, in localised tumour regions. Tumour cells exposed to such fluctuating hypoxia dynamics experience a selective advantage for malignant growth (Höckel et al. 1996) and resistance to chemo-and radiotherapy (Harrison and Blackwell 2004;Gray et al. 1953;Horsman et al. 2012). In spite of the obvious impact on tumour behaviour, the mechanisms and structural irregularities that contribute to oscillatory tumour blood flow remain unclear. We postulate that the existence of many redundant vessels in tumour networks, combined with the intrinsic nonlinearities of microscale blood flow, can play a significant role in such tumour blood flow fluctuations.
As a first step towards better understanding the microscale mechanisms that lead to unsteady flows and cycling hypoxia in tumours, we revisit the simplest model for network redundancy-a three-node network in which the two main flow branches are connected by a redundant vessel. Motivated by irregular tumour networks, where regulatory angiogenic mechanisms may be disrupted, and a range of vascular diameters may prevail, we consider different branch diameters in our model. We will show that such differences can have a significant affect on the existence of oscillatory dynamics. We combine numerical simulations and stability analysis to study the blood flow dynamics of the three-node network. By varying the branch diameter ratio (representing structural effects) and the inlet haematocrit level (representing the effect of local conditions), we construct a two-parameter stability diagram that delineates regimes in which multiple equilibria and oscillatory solutions exist. Using the haematocrit splitting model of Pries et al. (1989), we demonstrate how features of the model, particularly its nonsmoothness, affect the emergence of oscillatory instability.
In Sect. 2, we describe the model for blood flow in a three-node network and introduce the method we use to simulate its time evolution. In Sect. 3, we use dynamic simulations and linear stability analysis to characterise the steady and dynamic behaviours of the flow in the network. We summarise our conclusions in Sect. 4. Technical details relating to the analysis carried out in Sect. 3 are presented in "Appendix".

Model for Unsteady Blood Flow in a Three-Node Capillary Network
We study the unsteady flow of blood through a series of cylindrical capillaries that form a three-node network. The dependent variables are the nodal pressures, vessel flow rates, and haematocrit distributions. We start in Sect. 2.1 by formulating the coupled equations describing the flow dynamics in a single vessel as a function of its inlet and outlet conditions. Then, in Sect. 2.2, we formulate the three-node network model, describing how the blood flow and haematocrit in the different vessels are related at vascular junction points (internal nodes). A computational algorithm for simulating the time evolution of the flow in the network is given in Sect. 2.3.

Blood Flow Dynamics in a Single Vessel
We start by describing the flow in a single vessel. Since we focus on blood flow in microcapillaries, we neglect inertial effects and consider viscous flow in a cylindrical vessel whose length, L, is much larger than its diameter, D. Under the assumption of radial symmetry, the flow is assumed to follow Poiseuille's law which, after averaging over the vessel length, L, takes the form: where Q(t) is the total volumetric flow rate and is the pressure drop along the vessel (x 0 and x L denote the vessel entrance and end points, respectively; for straight vessels: L = x L − x 0 ). Additionally, represents the vessel-averaged haemodynamic resistance, and μ(t) denotes the vesselaveraged viscosity which is given by Following Pries et al. (1994), we assume that the apparent blood viscosity, μ, depends on the vessel diameter, D (given in units of microns) and the discharge haematocrit, H = H (x, t), represents the ratio of RBC flux to total flow rate, as follows: In Eq. (3), In order to calculate the average viscosity in Eq.
(2), the spatio-temporal distribution of haematocrit has to be evaluated. When considering haematocrit transport, we assume the case of a plug flow, i.e. the radially averaged velocity of RBCs is equal to the radially averaged plasma velocity; we assume that we can neglect the Fåhraeus effect (Fåhraeus 1929). While this assumption is made to simplify the analysis, we note that Karst et al. (2015) have shown that including the Fåhraeus effect does not significantly impact the system dynamics. Under these assumptions, the haematocrit distribution within each vessel is governed by the following one-dimensional advection equation: where the radially averaged blood velocity U (t) is given by By integrating Eq. (4) along the length of the vessel, it is straightforward to show that the vessel-averaged haematocrit, Similarly, averaging Eq. (3) and differentiating with respect to time, we deduce that the average viscosity, μ(t), evolves as follows: Equations (5) and (6)  Integrating Eq. (7) with respect to t, we have that where τ (t) is the time taken for haematocrit to propagate along a vessel of length L = x L − x 0 . The inlet haematocrit H (x 0 , t) depends on the haematocrit in the parent vessel(s) (see Sect. 2.2). Since the haematocrit is constant along the characteristic curves, if the inlet haematocrit in a specific vessel is known, then the corresponding outlet haematocrit is given by where τ (t) is defined implicitly by Eq. (8). In cases for which t 0 U (s)ds < L (the initial inlet haematocrit, H (x 0 , 0), has yet to propagate along the length of the vessel), we assign initial conditions for the outlet haematocrit, We note that if the velocity changes sign during a simulation, additional complexities can arise. For example, in some cases, the haematocrit may not have reached the vessel end-point before it starts to propagate backwards. These complexities will not be discussed here, for the sake of brevity, and because no such sign changes occurred for any of the simulations of the three-node network reported in the current work.

Three-Node Network Model
We consider unsteady blood flow within a three-node microcapillary network in order to derive a simple model to investigate how structural features of the network may invoke multistability and oscillatory dynamics. The three-node network consists of six vessels (numbered 1-6), all of length L, and has two inlets and a single outlet (see Fig. 1). Motivated by the irregular structure of tumour vasculature, we focus on the effect of asymmetry in the diameters of the two main flow branches. We assume that the diameters of vessels 2, 3, 5 and 6 are identical and denote by D their diameters. We assume further that the diameters of vessels 1 and 4 are identical and denote by α the ratio of their diameters to those of vessels 2, 3, 5 and 6. As we show below, this difference in vessel diameters plays an important role in the emergence of oscillatory dynamics in the three-node network.
Blood flow dynamics are greatly affected by the local haematocrit level due to its affect on blood rheology. While macroscale haematocrit levels are generally uniform in vascular networks (H ≈ 0.45 in humans), microscale haematocrit levels can be extremely heterogeneous due to plasma skimming effects (Pries et al. 1992). In order to study the combined effect of geometrical features and local haematocrit conditions, we assign a constant haematocrit discharge H 0 at both inlets. In this study, the inlet haematocrit, H 0 , and the diameter ratio, α, serve as the two key parameters governing the system dynamics.
We impose constant and equal pressure differences between both inlets and the outlet nodes: a discussion of alternative boundary conditions is included in Sect. 3.2.
In Sect. 2.1, we considered flow in a single vessel, assuming that the boundary conditions (pressure difference and inlet haematocrit) for this vessel could be determined from the flow in other vessels. We now formulate a system of algebraic equations for the nodal pressures based on mass conservation applied at internal nodes. Each internal node (marked with red circles in Fig. 1) represents a junction between three vessel segments (see Fig. 2). Accordingly, application of conservation of mass at an internal node yields where the indices a, b, and c refer to fluxes from the central node, where the pressure is p 0 , towards connected nodes with pressures p a , p b , and p c , respectively, as shown in Fig. 2. This notation is generic, i.e. it does not relate to a specific junction; later, we will use the notation Q i (i = 1, 2.., 6) to denote fluxes in specific vessels within the three-node network. Combining Eqs.
(1) and (10), we arrive at the following equation for the nodal pressure at a generic junction where represents the vascular conductivity [1/R j in Eq. (1)] as a function of the specific diameter, D j , and average viscosity, μ j , in the vessel connecting node 0 and the nodes marked with j = a, b, c. If we apply Eq. (11) at all internal nodes of the network, and prescribe the pressures of the inlet and outlet nodes, then we obtain a system of algebraic equations for the nodal pressures. In order to close these equations, knowledge of the average viscosity within all network vessels, μ i (i = 1, 2, . . . , 6), is required. In order to determine the viscosity, we must also calculate the corresponding discharge haematocrit values, H i (x, t). In Sect. 2.1, we studied haematocrit propagation through a single vessel. For a network, we must relate the haematocrit at the entrance to that vessel and the haematocrit at the end of the parent vessel (or vessels) that supply it with haematocrit. For example, consider the node connecting vessels 2, 3, and 5 in Fig. 1 (denoted as node 2-3-5 from now on). If the flows through vessels 2 and 3 converge into vessel 5, then the haematocrit at the entrance to vessel 5 is determined by applying a mass balance to the RBCs, Here, the vessel edges x 0 = 0 and x L = L refer to points in the coordinate system for a specific vessel ordered so that the pressure at x 0 is higher than that at x L . Alternatively, if the flux in vessel 3 changes direction, such that the flows in vessels 3 and 5 diverge from vessel 2, then a haematocrit splitting rule is applied, where the functions F 3 and F 5 = 1 − F 3 (due to haematocrit conservation) define the splitting rules for vessels 3 and 5, respectively. In the simplest nonlinear models, e.g. the model of Klitzman and Johnson (1982), the splitting rule depends solely on the flux ratio between the daughter and parent branches. More generally, however, the splitting rule may also depend on the branch diameters and haematocrit. In this work, we use the model of Pries et al. (1989). For the three-node network under consideration, haematocrit splitting only occurs at either node 2-3-5 or node 1-3-4, and we denote the splitting function at node 2-3-5 as F 3 = F(ψ), where ψ = Q 3 /Q 2 , while the splitting function at node 1-3-4 is F 3 = F * (ψ * ), where ψ * = Q 3 /Q 1 . For illustrative purposes, here we define the splitting function at node 2-3-5: where In Eq. (14), we have used the fact that the haematocrit in vessel 2 is equal to the inlet haematocrit, H 0 . The splitting function in node 1-3-4, F * (ψ * ), can be readily obtained by replacing D 2 by D 1 and D 5 by D 4 in the splitting function's coefficients A, B, and ψ 0 , following Eq. (14).
It is important to note that at bifurcations, the conserved quantities are the total flow, 3 i=1 Q i = 0, and the haematocrit flow, 3 i=1 Q i H i = 0, but not the haematocrit concentration. Due to the nonlinear form of the splitting rule, the proportion of haematocrit flux bifurcating into the favoured daughter branch may be larger than the proportion of total haematocrit and plasma flux bifurcating into this branch ('plasma skimming'), leading to increased haematocrit concentration in the favoured daughter branch relative to the parent branch. The nonlinearity of the splitting rule facilitates non-uniform distribution of haematocrit between the different vessels.
With all the components for modelling blood flow within a microcapillary network defined, we now explain how we construct model solutions for the nodal pressures, vessel flow-rates, and average haematocrit.
With the haematocrit distributions known at an initial time, the nodal pressures can be found using Eq. (11), and the flux (and average velocity) in each vessel can be determined using Eq. (1). Then, the inlet haematocrit for each vessel can be deduced from the outlet haematocrit of its parent vessel(s) (13) or (12). Evaluating the haematocrit propagation time, τ , using Eq. (8), and assigning to Eq. (9), the right-hand sides of Eqs. (5) and (6) are determined, such that the time evolution of the average haematocrit H and viscosity μ can be calculated. The numerical algorithm for implementing the outlined model is given below.

Dynamic Simulation Algorithm
The following algorithm was used to generate numerical solutions for the timedependent flow and haematocrit propagation in the capillary network: 1. The initial inlet haematocrit H i (x 0 , 0), outlet haematocrit H i (x L , 0), and average haematocrit H i (0), are prescribed for each vessel (i = 1, 2, . . . , 6), such that they are equal at a specific vessel. The corresponding initial values of the average viscosity, μ i (0), are then calculated in each vessel. A constant haematocrit value, H 0 , is prescribed in the inlet vessels. In such vessels: Constant pressures, p in = P and p out = 0, are prescribed at the inlet and outlet nodes, respectively, such that a constant overall pressure-difference, P, is imposed between both inlet nodes and the outlet node. 2. We adopt the methodology employed by Davis and Pozrikidis (2011) to solve Eq. (11) for the internal nodal pressures at each time step using a Gauss-Seidel iterative method. Iterations are continued until the relative change in value of each nodal pressure on subsequent iterations is less than a prescribed tolerance of 10 −20 . 3. The fluxes Q i and average velocities U i are calculated for all vessels (i = 1, 2, . . . , 6) using Eq. (1). 4. At each timestep (n = 1, 2, . . . , N , t n = n t; t denotes the magnitude of the time increment), and in each vessel (i = 1, 2, . . . , 6), the haematocrit propagation time, τ n i , is calculated by numerically integrating Eq. (8). The integral in Eq. (8) is evaluated for t − τ set equal to each discrete time in the interval 0, t n . The approximated value of t −τ is then determined as the time at which the left-hand side of Eq. (8) changes sign. When the flow in the vessel does not change its direction, the search for the value of t − τ can be restricted to the interval (t − τ ) n−1 where (t − τ ) n−1 i is the retarded time from the previous timestep. This feature is crucial in order to obtain reasonable run times of the simulation code. 5. The haematocrit at each vessel end point, H i (x L , t), is evaluated using Eq. (9). 6. Depending on the flow directions in the vessels connected to each internal node, either a haematocrit flux balance Eq. (12) or a splitting rule Eq. (13) is used to update the inlet haematocrit value, H i (x 0 , t), in the daughter branch(es). 7. The average haematocrit and viscosity values in each vessel are numerically advanced in time by applying Euler quadrature to Eqs. (5) and (6). 8. Steps 2-6 are repeated until t = t N .
The above algorithm was implemented in MATLAB. The code is available at the following GitHub repository: https://github.com/yaronbenami/blood_flow.

Characterising the Steady and Dynamic Behaviours of a Three-Node Network
In order to streamline the analysis of the three-node network considered in Fig. 1, we nondimensionalise the governing equations. Fluxes are normalised by the steady-state flow rate in vessel 2 so that where the superscript (0) denotes a steady-state value and, hereafter, tildes denote dimensional quantities. When pressure boundary conditions are prescribed, Q (0) 2 changes as the system parameters α and H 0 vary. However, this scaling was chosen because it renders the dimensionless formulation of the haematocrit splitting rule less cumbersome. Average vessel resistances, as defined by Eq. (1), are scaled by the constant resistance in vessel 2, where D = D/1 µm as appropriate for using the viscosity function of Pries et al. (1994) [Eq. (3)], α i = α for i = 1, 4 and α i = 1 otherwise (see Fig. 1). The spatial coordinate is scaled by L, such that for each vessel x ∈ [0, 1]; time is nondimensionalised by the time for haematocrit to propagate through vessel 2 when the flow is steady Under these scalings, the dimensionless parameters governing the system behaviour are α, H 0 , and D, and the parameter values used in this work are such that α ∈ [0.25, 2.25], H 0 ∈ [0, 1], and D = 20.

Steady-State Solutions
For the three-node network, if the same pressure is imposed at both inlets, then a steady solution with no flow in vessel 3 always exists. In this case, the haematocrit in all vessels is equal to H 0 , except for vessel 3, which has no haematocrit. It is straightforward to show that the dimensionless flux in the upper branch (scaled by 2 ) is given by Henceforth, we refer to this steady solution as the "trivial solution", and vessel 3 as the "redundant vessel", because there is always a steady state for which it transports neither haematocrit nor plasma. With Q In previous work, Gardner et al. (2010) showed that the three-node network admits multiple steady-state solutions. They also showed that the network possesses three solutions if the inlet haematocrit, H 0 , exceeds a threshold value. However, they did not characterise the stability of the steady-state solutions. Therefore, in this section, we characterise the multiple steady-state solutions of the network in order to subsequently study their stability characteristics (in Sect. 3.2). We start by formulating the equations that define the two non-trivial steady-state solutions: • Case I-flux flows from the bottom to the top branch (node 2-3-5 to node 1-3-4).
These solutions differ in the node at which the haematocrit splitting rule is imposed: for Case I, the haematocrit splitting rule is imposed at node 2-3-5, while for Case II it is imposed at node 1-3-4.
At a steady state, the haematocrit in each vessel is independent of spatial position x [we set ∂/∂t = 0 in Eq. (4)], and, thus, the resistance is also independent of x. Further, the pressure difference in vessel i depends on the flux via where from Eq. (1), If we consider the steady-state pressure drop along the loop formed by the three internal nodes, then we have that For Case I (i.e. flux flows from the bottom to the top branch), substituting (16) and (17) into (18) yields In Eq. (19), we have also exploited the mass balance Eq. (10) at nodes 1-3-4 and 2-3-5, which yields Q 3 , respectively. Additionally, imposition of equal inlet pressures yields Substituting (16) into (20), together with the prescription of inlet haematocrit H 0 , we find that Equations (19) and (21) should be supplemented by expressions for the haematocrit in vessels 3, 4 and 5. The haematocrit splitting rule (14) at node 2-3-5 can be written as Equation (22), together with haematocrit mass balance (12), yields the following expression for the haematocrit in vessel 5, while the haematocrit mass balance at node 1-3-4 yields the following expression for the haematocrit in vessel 4, Equations (19)-(24) define the nontrivial steady-state solution for Case I, where blood flows from the bottom to the top branch. For Case II, where blood flows from the top to the bottom branch, a similar analysis leads to the following system of equations: Here, we have exploited conservation of mass at nodes 1-3-4 and 2-3-5, which now means that We note that for both Cases I and II we assume Q (0) 3 > 0, thus the change of flux direction manifests via sign changes between Eqs. (19), (21) and (25), (26), respectively.
Equations (19)-(24) or (25)-(27) can be solved numerically to obtain the steadystate solutions in terms of the system parameters α and H 0 . Naturally, both sets of equations reduce to the trivial solution when Q (0) 3 = 0. As the trivial solution exists for all parameter values, a bifurcation should occur when one of the nontrivial solutions for Q (0) 3 approaches Q (0) 3 = 0. The strategy we use to identify such bifurcation points is to first find nontrivial steady-state solutions to Eqs. (19)-(24) and (25)-(27) for sufficiently large values of H 0 -one for each direction of flow in the redundant vessel. Then, we use numerical continuation to track these solution-branches as H 0 decreases. In Fig. 3, we present the multiple steady-state solutions found using this numerical tracking technique for α = 0.45. For all calculations, we assumed a nominal, dimensionless diameter of D = 20 (which corresponds to a dimensional diameter of 20 µm). As mentioned above, for both Cases I and II, Q 3 is considered positive.
3 > 0; red curve in Fig. 3), then the nontrivial solution branch ceases to exist at a saddle-node (fold) bifurcation at which H 0 = H S < H T (H S = 0.3288 in Fig. 3).
The thin dash-dotted black lines in Fig. 3 denote the critical fluxes, |Q 3 = 0 in the grey region bounded by these two lines. Following Karst et al. (2015), we term these critical values "skimming thresholds".
The function used in Eq. (14) imposes nonsmoothness of the splitting rule at the skimming threshold. The emergence of a fold bifurcation at one of the skimming thresholds (dash-dotted black lines in Fig. 3) might lead to the conjecture that the nonsmoothness of the haematocrit splitting function is responsible for the emergence of the fold bifurcation. However, a separate analysis, using a smooth splitting rule (see Sect. 3.4), yields a qualitatively similar bifurcation diagram which, for the sake of brevity, is not presented here. We conclude that the bifurcation structure shown in Fig. 3 is due to the network configuration, rather than the specific splitting rule.

Stability of the Trivial Solution
The bifurcation diagram presented in Sect. 3.1 suggests that the trivial solution exchanges stability with nontrivial steady-state solutions at a transcritical bifurcation point. To verify this finding, we conducted dynamic flow simulations, using the trivial solution to initialise the network flow. Figure 4 shows how the vessel-averaged haematocrit (Fig. 4a, b) and fluid flux in the redundant vessel (Fig. 4c, d) change over time when we fix the ratio of the vessel diameters and inlet haematocrit so that α = 0.5 and H 0 = 0.45. The simulations evolve to either a different (nontrivial) steady state (Fig. 4a, c) or an oscillatory state (Fig. 4b, d). The two different dynamics were obtained by imposing small perturbations (±10 −6 ) on the inlet haematocrit of vessel 2 at the first time step; a positive perturbation resulted in evolution to a nontrivial steady state (Fig. 4a, c), while a negative perturbation resulted in the onset of oscillatory dynamics (Fig. 4b, d). The insets in Fig. 4c, d illustrate the different flow directions in vessel 3. A positive flux in the redundant vessel Q 3 > 0, i.e. blood flows from the higher (α < 1) to the lower resistance branch, leads to a nontrivial steady state (corresponding to the red line in Fig. 3), while Q 3 < 0 may lead to oscillatory dynamics.
In light of these findings, we performed a linear stability analysis of the trivial steady-state solution in order to characterise its local stability. We made the following ansatz for the flux, average resistance, and haematocrit in the i-th vessel: where i , and H  We note here that the above perturbation equations are for Case I (blood flows from the bottom to the top branch, i.e. Q (0) 3 < 0 in Fig. 4). While the dynamic simulations have shown that the flow direction in the redundant vessel affects the attractor to which the system evolves, we will now show that it does not affect the stability of the trivial solution (the effect of the flow direction in the redundant vessel will be clarified in Sect. 3.3, where the stability of the nontrivial steady-state solutions will be examined).
To obtain the haematocrit perturbation in vessel 3, we state the haematocrit splitting rule at node 2-3-5: Substituting from (28) in (29), with Q (0) 3 = 0, and equating terms of O( ), we deduce that This is because when we use the model of Pries et al. (1989), F(Q (14)]. We note further that the trivial steady-state is only possible for haematocrit splitting models with F(Q (0) 3 = 0) and ∂ F/∂ψ| Q (0) 3 =0 = 0. Additionally, the linear stability of the trivial solution depends on the value of ∂ 2 F/∂ψ 2 | Q (0) 3 =0 , which is identically zero in Eq. (14). These features will be important when we discuss the effect of smoothing the splitting function in Sect. 3.4.
Applying the haematocrit mass balance, Eq. (12), at nodes 2-3-5 and 1-3-4, we find that In Eq. (31) we have used Eq. (15) from which we have that Q To complete the haematocrit distribution, we apply mass balance at node 4-5-6 to obtain Having defined the haematocrit perturbations in terms of the flux perturbation q 3 , we now relate the hydrodynamic-resistance perturbations to the haematocrit perturbations. In the unsteady case, the pressure difference is related to the flux via the average resistance [Eq. (1)], where the nondimensional average resistance reads, and Expanding R i as a regular power series in the small parameter 1, we obtain Eq. (17) at leading order, and at O( ), we have Then, we relate the fluxes and resistances in the different vessels by applying Eq.
(1) to the three constraints on the pressure drops in the network. At O( ), Eq. (18) reads where we have used the result from Eqs. (30) and (34) that r 3 = 0. Since we impose a constant pressure difference P between the inlet and the outlet nodes, we may assume, without loss of generality, that the O( ) perturbation to P is zero. Then, we have that where p ( ) i (i = 1, 2, . . . , 6) denote O( ) pressure differences. Substituting for p ( ) i from Eqs. (1) and (28) into Eq. (36), we deduce that and where, since we impose constant inlet haematocrit values (H 1 = H 2 = H 0 ), we have assumed, without loss of generality, that h 1 = h 2 = 0; this leads, via Eq. (34), to r 1 = r 2 = 0. Three additional equations are obtained by balancing the flow at each node: Equations (30)-(39) form a transcendental eigenvalue problem for λ. In practice, however, when perturbing about the trivial steady-state solution, λ attains only real values. In order to show this, it is helpful to consider the simpler problem of fixedflux boundary conditions. In this case, we fix q 1 = q 2 = 0 instead of imposing zero O( ) pressure drops between the inlet and outlet nodes [Eqs. (37) and (38)]. Then, the eigenvalue problem reduces to a single equation, where λ α = α 2 R 1 λ. For any choice of λ, the imaginary parts of the two terms on the right-hand side have the form where λ = σ + iω, θ = tan −1 (ω/σ ), and C and C α are positive coefficients. It can be readily shown that these two terms are both positive (when ω < 0), both negative (when ω > 0), or both zero (when ω = 0). Since the left-hand side of Eq. (40) has no imaginary part, it follows that the equation is only satisfied for real λ.
The same argument can be applied (although the calculations are more cumbersome) when pressure boundary conditions are imposed (results not shown). We conclude that the trivial solution cannot undergo a Hopf bifurcation and that transitions to the two attractors shown in Fig. 4 occur when the trivial solution exchanges stability with one of the nontrivial steady-state solutions, as shown in Fig. 3. With ω = 0, linear stability analysis of the trivial solution cannot explain the transition to the oscillatory state shown in Fig. 4b, d. In Sect. 3.3, we demonstrate that the onset of oscillatory dynamics occurs via a Hopf bifurcation from one of the nontrivial steady-state solutions. In order to find the critical conditions for instability of the trivial steady-state solution, we analyse Eq. (40) in the limit as λ → 0: With R The mechanisms driving instability of the trivial steady-state solution can be explained as follows. Suppose, without loss of generality, that a small flow perturbation, with zero haematocrit, enters redundant vessel 3 from junction 2-3-5 (as mentioned above, since the linear stability of the trivial solution does not depend on the flow direction in the redundant vessel, it suffices to consider this case). Consequently, the haematocrit and resistance in vessel 4 decrease while those in vessel 5 increase. This leads, at supercritical conditions, to more flow being redirected towards vessel 4, creating a positive feedback mechanism which destabilises the trivial solution.
At the critical conditions (λ = 0) given by Eq. (41) (when flux boundary conditions are imposed), the perturbations in the pressure-drop along the three internal nodes (vessels) caused by the increase in the resistance of vessel 5 and the corresponding decrease in the resistance of vessel 4 [the terms in brackets in Eq. (41)] are balanced by the pressure drop due to the flow perturbation in vessels 3, 4, and 5. This condition yields a relationship between a particular value of α and the critical H 0 for instability of the trivial solution. For inlet haematocrits that are larger than the critical value, this balance cannot hold due to the increase in the resistance perturbation, resulting in the positive feedback mechanism described above.
We sketched the critical curves by discretising α in the range [0.25, 2.25] and solving numerically for H 0 (α) using a continuation scheme, where α is a continuation parameter. In Fig. 5, we plot the critical curves thus obtained for the two types of boundary conditions. To confirm that the critical curves correspond to the transcritical bifurcation points, H T , evaluated in Sect. 3.1, we calculated discrete values of H T (α) by locating the bifurcations in the steady-state solution diagrams for a range of values of α. Comparison with the critical curve in Fig. 5 indicates good agreement, providing independent validation of our stability analysis. Figure 5 shows that for a given value of α, the critical haematocrit for the fixedpressure condition is always larger than the corresponding value for constant-flux conditions, suggesting enhanced instability of the latter. From a mechanistic perspective, we note that under the assumption of Poiseuille flow, the flow rate is proportional to the pressure gradient, This means that imposing a flux condition is equivalent to applying a condition on the derivative of the pressure. Typically, a condition on a derivative is less restrictive than a condition on the variable itself, leading to a reduced parameter range of stability. Our findings are consistent with experimental results reported by Storey et al. (2015). They conducted experiments involving mixtures of two fluids in slightly different network geometries and found that imposing flux boundary conditions produced a larger region of instability than imposing pressure boundary conditions. We now analyse the stability diagram in the limits of α 1 (α 1), representing the cases when the diameter of the top branch is much smaller (larger) than the diameters of the bottom branch and the redundant vessel.

Limit of˛ 1
From Fig. 5, we note that the critical stability curves for both types of boundary conditions are similar when α 1. In this limit, R Additionally, R When flux boundary conditions are imposed, Eq. (41) reduces to When pressure boundary conditions are imposed, analysis of Eqs. (30)-(39) in the limit α 1 shows that a nonzero solution is only possible if the flux perturbations satisfy the following scaling Assigning scaling (46) to equations (30)-(39) in the limit α 1 yields an equation which is identical, at leading order, to Eq. (45). The dashed black line in Fig. 5 corresponds to the solution of Eq. (45); this line is indistinguishable from the two solid curves for α 0.4.

Limit of˛ 1
It is clear from Fig. 5 that for both types of boundary conditions, the critical curves tend to different constant values of H 0 when α 1. In this limit, the steady-state resistance in the branch with diameter α D is small: and the perturbation to the resistance in this branch is much smaller than the perturbation in the branch with diameter D: Therefore, all the terms which include α in Eq. (41) for fixed-flux conditions (or Eqs. (30)-(39) for fixed-pressure conditions) are negligible when α 1, rendering the two sets of equations independent of α, a result which is consistent with the numerical results in Fig. 5. When flux boundary conditions are imposed, Eq. (41) reduces to When pressure boundary conditions are imposed and α 1, analysis of Eqs. (30)-(39) shows that a nonzero solution is only possible when the flux perturbations are of similar magnitude, that is, Assigning the scaling in Eq. (50) to Eqs. (30)-(39) yields The black dash-dotted lines in Fig. 5 correspond to the asymptotic values of the critical haematocrit, H 0 , when α 1. For both types of boundary conditions there is excellent agreement between the asymptotic values of the critical inlet haematocrit for instability and the values given by Eqs. (41) or (30)-(39). We note further that since the left-hand sides of Eqs. (49) and (51) are monotonically increasing in H 0 , the asymptote for the pressure boundary conditions is larger than that for the flux boundary conditions.

Oscillatory Solutions
Given the oscillatory dynamics predicted by the dynamic simulations in some cases, it is of interest to identify regions of parameter space in which they exist. To this end, we performed linear stability analysis of the nontrivial steady state solutions where, as hinted by Fig. 4, the flux in the redundant vessel flows from the lower to the higher resistance branch (this assumption will be justified a posteriori). In order to cover the full range of solutions, we consider separately Cases I and II from Sect. 3.1. As demonstrated in Sect. 3.1, the key difference between these cases is the node at which the haematocrit splitting rule is imposed.
We use the same ansatz as in Eq. (28) to perturb the steady-state equations. While the full analysis is presented in "Appendices A and B", important differences between perturbations to the trivial and nontrivial states are emphasised here. For Case I, when blood in the redundant vessel flows from the bottom to the top branch (analysis for Case II is presented in "Appendix B"), Eqs. (19)-(24) define the steady-state solutions. Linearising the splitting function at node 2-3-5 [Eq. (29)] about its steady state [Eq. (22)], we deduce that In contrast to the perturbation of the trivial state, here Q (0) 3 = 0, leading to h 3 = 0. The haematocrit mass balance at node 1-3-4 yields the following expression for the perturbation to the haematocrit in vessel 4: (53) In Eq. (53), we note the exponential term on the right-hand side, which represents the effect of the time-delay between haematocrit entering vessel 3 and its propagation into vessel 4, induced by the nonzero value of Q (0) 3 . This time delay plays a key role in the emergence of oscillatory dynamics.
The equations governing the perturbations v = {h 3 , . . . , h 6 , r 3 , . . . , r 6 , q 1 , . . . , q 6 } T can be written as a linear system of the form where the matrix A(λ, V (0) ) depends on the eigenvalue, λ = σ + iω, and the steadystate solution, 3 , H 3 , H 4 , H 5 }, which satisfies Eqs. (19)-(24) for Case I (or Eqs. (25)-(27) for Case II). Full statements of the equations governing the perturbations v are presented in "Appendices A" (Case I) and B (Case II). Equation (54) constitutes a transcendental eigenvalue problem. In practice, we determine the system's iso-σ (iso-growth-rate) contours in the (α, H 0 ) parameter space by seeking solutions for {H 0 , ω} as a function of α that satisfy for specific values of σ . To carry out the bifurcation analysis, we used a numerical continuation scheme to obtain iso-σ curves in the (α, H 0 ) parameter space, with α as the continuation parameter. We initialised the continuation scheme using solutions estimated from the dynamic simulations presented in Sect. 2.3 (for a given combination of values of α and H 0 , we estimated ω and σ by evaluating the oscillation frequency and growth rate within a short time period after the bifurcation commenced). The resulting stability diagram is presented in Fig. 6, for iso-σ curves corresponding to Cases I (α < 1) and II (α > 1), respectively. The solid black line marks the critical curve for stability of the trivial steady-state solution and the dashed black line represents the skimming threshold of the redundant vessel. In the grey region that separates the critical stability curve of the trivial steady-state solution and the skimming threshold, there is no haematocrit in the redundant vessel, although |Q 3 | > 0. Remarkably, in this case, the nontrivial steady states are stable; none of the contours with positive growth-rates cross the skimming threshold, indicating that oscillatory instability can only occur for values of inlet haematocrit H 0 above this threshold. The requirement for haematocrit to be present in the redundant vessel in order to generate oscillatory dynamics can be attributed to the effect of time-delays in the system: when no haematocrit is present in the redundant vessel (H (0) 3 = h 3 = 0), perturbations in Q 3 lead to instantaneous changes in the haematocrit (and, consequently, the resistance) in vessel 4 (or vessel 5 when the flow in vessel 3 is in the opposite direction). Thus, we conjecture that in the absence of time delays, self-induced oscillations cannot be sustained. Figure 6 suggests that oscillations only occur in the presence of multiple equilibria. (The oscillatory regime is always above the critical curve of the trivial solution, indicating the presence of other nontrivial steady-state solutions.) However, the existence of multiple equilibria is not a necessary condition for the existence of oscillatory states. For example, Karst et al. (2015) studied a slightly different network geometry and identified small regions of parameter space in which oscillatory solutions exist in the presence of only a single steady-state solution. These regions, however, seem to exist only when the diameter of the redundant vessel is very small.
In Fig. 6, the iso-σ contours do not form closed contours, but they originate from the skimming threshold. We conclude that there is a jump between negative and positive growth rates at the skimming threshold. This discontinuity arises because there is a nonsmooth Hopf bifurcation-a consequence of the singularity in the system's Jacobian at the skimming threshold, which arises when the exponent B in Eq. (14) is such that B < 2. For the parameter regime used in Fig. 6, B ∈ [1.14, 1.37] on the skimming threshold and, hence, we have a nonsmooth bifurcation (the effect of smoothing the splitting function is discussed in Sect. 3.4). The nonsmooth Hopf bifurcation at the skimming threshold was reported by Karst et al. (2015), who determined the critical conditions for oscillations in a different network geometry, using a similar haematocrit splitting model.
Interestingly, Fig. 6 reveals an additional stable region in a neighbourhood of α = 1 (i.e.f where the vessel diameters in the two branches are similar), suggesting that a critical difference in the flux between the two branches is needed to trigger oscillatory solutions. Recall that the two regions in which oscillatory solutions exist in Fig. 6 correspond to the two flow directions in vessel 3 (oscillations in the bottom-to-top flow configuration are restricted to α 1, and to α 1 when the flow is in the opposite direction). The fact that the critical curves for the onset of oscillations (σ = 0) in Fig. 6 do not cross the line α = 1 demonstrates that, in this network geometry, oscillatory solutions can only exist when the flux in the redundant vessel goes from the lower to The oscillatory dynamics can be understood by considering the two sources of nonlinearity: (i) the haematocrit-dependent viscosity of blood (i.e. the Fåhraeus-Lindqvist effect), and (ii) the nonlinear splitting of haematocrit at vessel bifurcations ("plasma skimming"). While the former induces coupling of the flow and haematocrit concentration, the latter allows non-uniform haematocrit distributions throughout the network. Both effects are essential for the feedback mechanisms that generate self-sustained oscillations: plasma skimming leads to relatively little haematocrit entering the redundant vessel which, in turn, dilutes the haematocrit and, consequently, reduces the resistance to flow (due to the Fåhraeus-Lindqvist effect) in the smaller-diameter branch, so that more flow is redirected to the redundant vessel. The time-delayed negativefeedback is also a consequence of the haematocrit-dependent viscosity: the increase in haematocrit in the redundant vessel leads to a delayed increase in the resistance of this flow path. Since no other sources of nonlinearity are included in our model, we believe that it represents a minimal model of how the presence of redundant vessels can promote oscillatory blood flow when physiologically realistic rules are used to describe the Fåhraeus-Lindqvist effect and plasma skimming.
Consider, for example, the feedback between the flow and haematocrit in vessel 3 (the redundant vessel) and vessel 4 when the flow in vessel 3 is directed towards vessel 4 (oscillatory solutions for α 1 and steady solutions for α 1). The time evolution of the fluid flux and haematocrit in this case are presented in Fig. 7 for H 0 = 0.5 and two values of α < 1. The following positive feedback mechanism acts on the flux in vessel 3: when Q 3 is increasing, more haematocrit enters vessel 3 (if the state of the system is located above the skimming threshold), and less haematocrit enters vessel 4 (see the counter-phase behaviour of Q 3 and H 4 in Fig. 7). In effect, the increase in Q 3 dilutes the haematocrit in vessel 4, while there is a time delay before the increase in the haematocrit entering vessel 3 propagates through the vessel and reaches vessel 4 (see the phase-lag between Q 3 and H 3 in Fig. 7). The decrease in haematocrit reduces the resistance in vessel 4, so that eventually more flow is redirected from the bottom to the top branch, through vessel 3. This half cycle is reversed when the increase in the resistance of vessel 3 (due to the increase in its haematocrit) becomes large enough to reduce Q 3 (see correspondence between the maxima of H 3 and the times at which Q 3 starts to decrease in Fig. 7). While similar arguments can be applied when α > 1, changes in the haematocrit entering vessel 4 in response to perturbations in Q 3 are governed by the ratio of the steady-state fluxes Q  Fig. 7b, d) and, thereby, suppresses the positive feedback mechanism.
In this work, we chose the nominal diameter to be D = 20 (20 µm in dimensional units) as a representative diameter of physiological microcapillaries. Some quantitative differences should be expected when D is changed. The critical value of H 0 at which oscillations emerge is dominated by the skimming threshold. When the nominal diameter D is decreased, the dominant effect is a reduction in the value of α at which the critical H 0 is minimized (α ≈ 0.55 for D = 20). This happens because, for a given haematocrit concentration, the change in resistance with the change in diameter (∼ d μ(H , D)/D 4 /dD) increases as the diameter decreases. Therefore, the diameter difference between the two branches (i.e. the degree of structural asymmetry) decreases. By contrast, the size of the stable region in a neighbourhood of α = 1 increases as D decreases. In this case, the resistance of the redundant vessel (diameter D) increases and, therefore, a larger diameter ratio is needed to drive flow through the redundant vessel.
While all results presented in this work are for equal inlet pressures and haematocrits, our calculations (not presented here for brevity) suggest that oscillatory solutions, having qualitatively similar dynamics to the oscillations presented here, can still exist when there is a relative inlet-pressure and haematocrit difference of a few percent. This shows that the oscillatory instability is not a unique feature of threenode networks with identical inlet conditions; rather it exists for a range of boundary conditions. The susceptibility of the three-node network to oscillatory dynamics for non-equal boundary conditions should help to realize these oscillations experimentally, as identical inlet conditions are technically challenging to achieve in practice.

Oscillations with a Smooth Haematocrit Splitting Function
Motivated by the nonsmooth Hopf bifurcation to an oscillatory state encountered when we use the haematocrit splitting model of Pries et al. (1989), it is natural to ask how smoothing the splitting function will affect the stability of the system. Therefore, we use a third-order polynomial to smooth the discontinuities in Eq. (14) and consider an alternative splitting function of the form where y L (ψ) and y U (ψ) are third-order polynomials and ψ L and ψ U represent points at which these polynomials intersect the original splitting function. A third-order these conditions are imposed in order to preserve the stability properties of the trivial solution. We also require so that the governing equations are smooth throughout the parameter space. The method used to construct smoothing polynomials that satisfy the above conditions is described in "Appendix C". Figure 8 shows how the stability of the oscillatory solutions changes when the splitting function in Eq. (14) is replaced by Eq. (56). The results for the original, nonsmooth function are also included to facilitate comparison. Noticeably, the smooth splitting function produces closed iso-σ contours that differ from those obtained using the original, nonsmooth splitting function. Comparing the dashed and solid blue curves in Fig. 8, smoothing the splitting function appears to increase stability for α 1.3, while for larger values of α this trend is reversed. As expected, sufficiently far from the skimming threshold, the smooth and nonsmooth haemaotcrit splitting rules yield stability results that are practically indistinguishable.
Considering the perturbed equations of the nontrivial steady-state solutions, the haematocrit splitting rule appears only in Eq. (52), where it attains the following functional form Therefore, the magnitude of S(Q 3 ) determines the size of h 3 /q 3 , the ratio of the haematocrit to the fluid flux in the redundant vessel. As such, it plays a major role in the positive feedback mechanism (see discussion in Sect. 3.3), associated with the onset of self-sustained oscillations.
The inset in Fig. 8 shows how differences between the smooth and nonsmooth splitting functions affect the behaviour at the point (α, H 0 ) = (0.58, 0.391). This point is located just above the skimming threshold (unstable using the nonsmooth splitting rule), but below the curve on which σ = 0 for the smooth splitting rule (stable). The function S(ψ) shows how the choice of haematocrit splitting rule affects the linear stability of nontrivial steady-state solutions and, therefore, how smoothing the splitting rule may change the system dynamics. The very large gradient of S(ψ) in the nonsmooth model above the skimming-threshold also explains why the solution becomes unstable immediately above this critical value. The different values of the steady-state fluxes Q 3 ) than the system with the smooth splitting function, which results in an oscillatory instability of the former, while the latter is stable.
From the experimental viewpoint, the functional form of the haematocrit splitting rule for low flow rates (in a neighbourhood of the skimming threshold) is challenging to measure because of the large noise-to-signal ratio that prevails at low flow rates. The stability results presented here for the cases of smooth and nonsmooth haematocrit splitting rules delineate the "boundaries" of possible behaviours for the given network. We postulate that other functional forms for smoothing the model of Pries et al. (1989) that satisfy the conditions specified in (57) and (58) will yield bifurcation curves within these two "bounds".

Conclusion
In this paper, we studied microcapillary blood flow in a three-node network, exploring its multiple equilibria and the transition to oscillations via dynamic simulations and stability analysis. While multiple steady-state solutions and self-sustained oscillatory solutions in microcapillary blood flows have been reported previously (see, for exam-ple, Karst et al. 2015 and refs. cited therein), to our knowledge, the microstructural characteristics that promote unsteady behaviour have not previously been identified. In this work, we have demonstrated that specific structural abnormalities, in the form of redundant vessels which connect two flow paths with different resistances, are key to the emergence of oscillations. We have clarified the feedback mechanisms, arising due to the coupling of these structural features with the intrinsic nonlinearities of blood flow at the microscale (i.e. Fåhraeus-Lindqvist effect and plasma skimming), which gives rise to oscillatory dynamics.
In our analysis, we defined a vessel as redundant if there is an equilibrium solution having zero-flow through that vessel. Such a "trivial" solution is typically unstable, with nontrivial steady-state solutions bifurcating (stable at the bifurcation point) for sufficiently large inlet haematocrit values. Remarkably, as we demonstrated using dynamic simulations, starting from the trivial state, the system may evolve to either a different steady-state solution or an oscillatory solution. The paths leading to these long term solutions are sensitive to small changes in the inlet conditions, which dictate the direction of flow in the redundant vessel. The sensitivity of the system to small fluctuations in the boundary conditions may lead to highly unstable behaviour if such a motif is embedded in larger networks. Additionally, we postulate that the maximum number of possible steady-state solutions should rise dramatically as the number of redundant vessels in a network increases, because each redundant vessel may support three solutions (no flow and/or flow in either direction). The large number of equilibrium states, together with sensitivity to small fluctuations, may explain why highly irregular, almost chaotic flow is a characteristic feature of many vascular tumour networks (Kimura et al. 1996;Brurberg et al. 2007;Gillies et al. 2018).
To quantify the critical conditions for instability as the ratio of branch diameters (representing the structural driving force) and inlet haematocrit (representing the effect of local flow conditions) vary, we performed stability analysis of the trivial solution. We found that the transition from the trivial steady-state solution to oscillations occurs in two steps-the trivial state loses stability to a nontrivial steady-state solution which, in turn, undergoes a Hopf-bifurcation. By performing linear stability analysis of the nontrivial steady-state solutions, we showed further that the combined effects of a redundant vessel and vessels that offer different resistances to flow (via different diameters in this work) is key to the emergence of self-induced oscillations. Also, we identified a feedback mechanism that facilitates the onset of oscillations; here, the diameter ratio between the two branches (affecting the flux in the redundant vessel) and the presence of haematocrit in the redundant vessel (allowing for time-delay in the system) are crucial ingredients for such positive feedback to occur. In future work, we aim to evaluate the effect of redundant vessels in larger vascular networks and to explore different motifs which may generate larger feedback loops and, thereby, larger scale oscillatory dynamics. Such an investigation should consider the coupled behaviour of multiple sources of oscillations, and how their frequencies and amplitudes are modulated. Studying the haematocrit oscillations in larger networks will ultimately enable us to evaluate their effect on tissue oxygenation, which is of considerable importance in understanding the process of cycling hypoxia in tumours (as mentioned in Sect. 1).
Traditionally, studies of blood flow in large networks did not consider in detail what type of boundary conditions should be imposed because, in general, the appropriate choice of boundary conditions for a microcapillary network is unknown. Fry et al. (2012) showed that the choice of boundary conditions imposed on large-scale microcirculatory networks can significantly influence the steady state flow rates. While most of the analysis in this paper was performed for constant pressure boundary conditions, we showed that changing to fixed-flux boundary conditions can destabilise the system, by reducing the critical inlet haematocrit at which the trivial solution becomes unstable. Therefore, in future work, it would be of interest to examine how the stability of larger networks (where there are many more internal nodes than boundary nodes) is affected by changes in the type of boundary conditions imposed.
In this study, we used a haematocrit splitting rule due to Pries et al. (1989); this model includes a threshold value of the daughter-to-parent flux ratio, such that haematocrit only enters the daughter branch if the flow rate exceeds this critical value (the "skimming threshold"). The skimming threshold gives rise to a discontinuity in the splitting rule which, in turn, gives rise to nonsmooth stability diagrams. We used spline-smoothing to eliminate the points of discontinuity in the model. In so doing, we obtained solutions which exhibited a smooth transition in parameter space between steady and oscillatory states, while also converging to the results of the nonsmooth model sufficiently far from the skimming threshold. In contrast to the steady-state solutions, the Hopf-bifurcation patterns are sensitive to small changes in the haematocrit splitting rules in the regime when the daughter-to-parent flux ratio is small (in the neighbourhood of the skimming threshold). This sensitivity of the emergence of oscillatory dynamics to small changes in the haematocrit splitting rule used introduces significant challenges regarding how to measure and model the haematocrit splitting that occurs at such low flow rates. which allow us to calculate the polynomial coefficients as a L = F(ψ L ) ψ 3 L and a U = F(ψ U ) (ψ U − 1) 3 .