Multiscale modelling of blood flow in cerebral microcirculation: Details at capillary scale control accuracy at the level of the cortex

Aging or cerebral diseases may induce architectural modifications in human brain microvascular networks, such as capillary rarefaction. Such modifications limit blood and oxygen supply to the cortex, possibly resulting in energy failure and neuronal death. Modelling is key in understanding how these architectural modifications affect blood flow and mass transfers in such complex networks. However, the huge number of vessels in the human brain—tens of billions—prevents any modelling approach with an explicit architectural representation down to the scale of the capillaries. Here, we introduce a hybrid approach to model blood flow at larger scale in the brain microcirculation, based on its multiscale architecture. The capillary bed, which is a space-filling network, is treated as a porous medium and modelled using a homogenized continuum approach. The larger arteriolar and venular trees, which cannot be homogenized because of their fractal-like nature, are treated as a network of interconnected tubes with a detailed representation of their spatial organization. The main contribution of this work is to devise a proper coupling model at the interface between these two components. This model is based on analytical approximations of the pressure field that capture the strong pressure gradients building up in the capillaries connected to arterioles or venules. We evaluate the accuracy of this model for both very simple architectures with one arteriole and/or one venule and for more complex ones, with anatomically realistic tree-like vessels displaying a large number of coupling sites. We show that the hybrid model is very accurate in describing blood flow at large scales and further yields a significant computational gain by comparison with a classical network approach. It is therefore an important step towards large scale simulations of cerebral blood flow and lays the groundwork for introducing additional levels of complexity in the future.


Introduction
Blood flow through cerebral microcirculation is the primary driver of oxygen transport and waste removal in the cortex, making microcirculation fundamental to cerebral physiology [1][2][3]. Yet microvascular networks have been given much less emphasis than other fascinating components of the brain, such as the white matter and the activity of neural networks [4]. Thus, we comparatively understand little about the structure and functions of these networks, let alone their regulation mechanisms, robustness to perturbations and involvement into pathologies. For instance, obstructions of small vessels during strokes, long-term capillary rarefaction or lymphocyte stalling occurring early in Alzheimer's Disease (AD) [5] may limit the transport of oxygen to the cerebral tissue and result in energy failure and neuronal death. Understanding how these disease-induced modifications of the network architecture affect transport mechanisms in the brain is therefore a fundamental challenge that may ultimately lead to clinical breakthroughs. One of the most complex aspect of the problem is that it involves a broad range of spatial scales [6], ranging from the scale of the whole brain (*10 cm, see ). This multiscale architecture makes it difficult not only to understand the mechanisms occurring at the different levels, but also to understand and assess the systemic impact of localized effects.
Advances in both imaging techniques and modelling approaches have made progress possible in the last decade. In particular, multiphoton laser scanning microscopy [10,11], a cuttingedge fluorescence microscopy with unprecedented microscopic scale resolution and cortical penetration depth, has enabled the spatio-temporal investigation of hemodynamics and mass transfers at macroscopic scale in the living brain of healthy or diseased rodents (e.g. [1,3,12,13]). In parallel, hemodynamically-based perfusion techniques are widely used in clinics to explore much larger volumes of human brain, typically the scale of the whole brain, with a resolution of approximately 10-100 mm 3 . With regard to imaging techniques, the investigation of μm-thick cortical section, where blood vessels have been injected with India ink for contrast enhancement [7]. (b) Macroscopic scale (*18 mm 2 × 300 μm): reconstruction of parts of the collateral sulcus by confocal laser microscopy [8]. (c) Mesoscopic scale (*5 mm 2 × 300 μm): region of interest in which vessels of more than 10 μm in diameter are colored in black and vessels of less than 10 μm in diameter are colored in red (diameters have been multiplied by 2 for visualization). In contrast with the capillary bed, the arteriolar and venular trees have a quasi-fractal structure [9]. (d) Microscopic scale (*0.07 mm 2 × 300 μm): detailed view of the capillary bed. The capillary bed is dense and space-filling over a cut-off length of * 50 μm [9]. the whole range of scales, from the microscopic to the brain scale, is therefore possible. For models, the most successful approaches are network representations that have been developed to simulate blood flow and mass transport in brain microcirculation [14][15][16][17]. However, contrary to imaging techniques, the flow models available do not cover the whole spectrum of spatial scales and simulations of transport mechanisms in the brain cannot yet be performed at clinically relevant scales.
To understand better the limitations of network models, we first need to describe them in more detail. Such approaches treat the vasculature as a network of interconnected tubes. The advantage of doing so is that the flow equations do not have to be computed over a complex three-dimensional geometry, as the relationship between average blood flow and pressure drop in each tube can be described analytically. This simplification of the geometry greatly reduces the computational cost, often yielding a simple linear system with relatively fewer variables compared to a complete three-dimensional resolution of the flow. The geometry and topology of the vascular network can be either extracted from available anatomical databases (1-10 mm 3 ) or synthetically generated to investigate larger volumes, as was proposed in [18,19]. Non-linearities induced by the complex rheology of blood in individual vessels (e.g. Fåhraeus, Fåhraeus-Lindqvist, and phase separation effects) can also be taken into account using empirical laws, and solved using iterative methods [20]. The main limitation of this framework is that the computational cost of the simulations quickly grows with the number of vessels. This is because, although it simplifies the geometry, the approach remains a microscopic scale description of the flow. Thus, simulations become intractable in large microvascular networks, especially when performing many statistical realizations (e.g. for studying the effect of capillary occlusions), when needing thousands of iterations to accurately capture nonlinear effects, or when coupling blood flow with oxygen transport. To the best of our knowledge, the most advanced computations have reached a volume of about 30 mm 3 (*250 000 vessels, [18]) so far, a number that is several orders of magnitude lower than what would be necessary for brain scale simulations.
Can we develop models with the potential to simulate transport at the scale of the whole brain? The simplest way to proceed would be to use homogenized descriptions whereby the velocity and pressure fields are averaged in space over a large number of vessels and the detail of the architecture at microscopic scale is not fully described. These microscopic details are implicitly represented via effective parameters (e.g. permeability) that characterize the flow at a larger scale, which we term mesoscopic scale (Fig 1(c)). Such approaches, which we term "homogenized" or "continuum approaches", are prominent in porous media sciences [21][22][23], the most famous example being Darcy's law. One of the advantages of such descriptions is that they lend themselves very well to the assimilation of experimental or numerical data obtained in relatively small volumes. For example, the permeability can be obtained by numerical simulations at microscopic scale using anatomical data acquired in only~1 mm 3 [14,24]. This class of approaches also significantly reduces the computational cost compared to network models, as computations are performed over coarse mesoscopic scale grids. Unfortunately, homogenized models are not applicable to the flow problem in the whole brain microvasculature. The primary reason is that simple average representations are only available when the structure exhibits a discrete hierarchy of scales, not a continuous one as the brain microvasculature (Fig 1). This continuous hierarchy of scales stems from the architecture of feeding arteriolar and draining venular trees, which exhibit a quasi-fractal geometry [9] (highlighted in black in Fig 1(c)).
To bridge the gap, new modelling paradigms are needed. One such paradigm is a hybrid approach where the capillary bed (highlighted in red in Fig 1(c) and 1(d)), which is homogeneous and space-filling above a cut-off length of *50 μm [9], is described using an averaged description, whereas feeding arteriolar and draining venular trees (highlighted in black in Fig 1(c)) are treated using the classical network approach. The term "hybrid" is a reference to the fact that this model is a hybrid between network and homogenized approaches. Such methods have already been introduced in the context of cardiac [25,26] or tumor [27] microcirculation. A major difficulty of such descriptions is to properly couple both frameworks, an issue that has not been addressed in these previous papers. The coupling is difficult primarily because variables that are seemingly identical do not really have the same physical meaning. For instance, the blood pressure in the continuum approach represents a field that is spatially averaged at mesoscopic scale, while the pressure of the network approach represents a crosssectional average within each tube. Simply imposing a continuity of the pressures at each coupling point, i.e., at each connection of an arteriolar or venular vessel to a capillary (red dots in Fig 2(a)), may therefore not accurately describe the strong gradients that occur locally and important errors could propagate through the whole system. Similar issues apply for concentrations or temperatures when molecular exchanges or heat transfers are considered.
In this paper, our goal is to devise a theory for coupling network and continuum descriptions of blood flow in the brain microcirculation, which lays the groundwork for modelling microvascular blood flow at the scale of the whole brain and extending to mass transfers. We focus here on blood flow with the following main simplifications: the capillary bed is assumed to be isotropic, homogeneous and its architecture is represented by simplified three-dimensional lattices, in the same spirit as in the two-dimensional study in [14]; further, the non- trees (black networks) plunge into the cortical tissue where their endpoints are connected to the capillary bed (light red). Because of the multiscale architecture of the human brain (Fig 1), we adopt two separate strategies to model blood flow at macroscopic scale, depending on the hierarchical position of the vessels in the microvascular network. A network approach is used in arteriolar and venular trees while the capillary bed is considered as a continuum. The continuum model is then discretized using coarse mesh cells (green) at mesoscopic scale, with size h x × h y × h z ' (250 μm) 3 . The side of a discretization cell is significantly larger than the capillary characteristic length l cap = 50 μm. The red dots correspond to the coupling points connecting arterioles (or venules) to capillary vessels, where coupling conditions must be applied. (b) We use anatomically accurate data to illustrate spatial correlations between the coupling sites. Minimum distances have been computed between all coupling points according to their types (arteriolar to arteriolar, venular to venular and arteriolar to venular). The histogram plot presents the percentages of couplings associated to the same range of minimum distance. linearities induced by red blood cell partitioning at microvascular bifurcations are neglected, following [1,14,16,18,28], so that the hematocrit is uniform. Based on these simplifications, we introduce an analytical model describing the coupling between the network (Section 2.2) and the continuum (Section 2.3) approaches. This model takes the form of a function linking the discrete pressures at the arteriolar and venular coupled endpoints to the cell pressures in a finite volume formulation of the homogenized equations. To derive this model (Section 2.4), we use a decomposition of the pressure field around each coupling point into two analytical parts: a linear variation that captures the local network structure of the underlying capillaries close to the coupling point and a spherical solution of the continuum equations further away. Our method is validated by comparing simulations obtained using this hybrid approach to a complete network description where all the vessels including the capillaries are described by the network approach. We first assess the robustness of the approach in various basic situations, with a small number (< 2) of coupling points (Section 3.1). We then go on to consider anatomically accurate arteriolar and venular trees that interact through many coupling points with an idealized capillary bed (Section 3.2). Besides a significant improvement in computational performance (Section 3.3), our results show that properly coupling the different modelling components in the hybrid approach-not just continuity of pressures at the coupling points-is crucial to derive accurate representations of pressure and velocity fields in the cerebral microcirculation. By thoroughly exploring the theoretical framework in model situations, we lay the groundwork for the development of more general models accounting for additional levels of complexity in the future.

Models and method
Here, we present the models describing mass and momentum transport in the microcirculation for the network approach (black arteriolar and venular trees in Fig 2), the continuum representation (coarse green mesh in Fig 2) that replaces the capillary bed (light red vessels in Fig 2) and the coupling conditions at the interface between these two frameworks. For consistency, the colors conventions introduced in

Main simplifications and hypotheses
The simplifications mentioned in the Introduction (isotropy and homogeneity of the capillary bed at mesoscopic scales, uniform hematocrit) are needed to introduce the analytical coupling model presented in detail in Section 2.4. In addition to these, we use several other hypotheses that are common to all frameworks. First, blood typically flows at about 10 mm Á s -1 in the arteriolar and venular trees and at about 1 mm Á s -1 in the capillary network, so that the Reynolds number is everywhere significantly lower than 1 [29] and inertial effects can be neglected. Second, the pulsatility is negligible because the Womersley number is lower than 0.1 [30]. We further assume that the mass density of blood ρ and the gravity vector g are constant in the brain microcirculation, so that the pressures used in all equations are defined as the fluid pressures minus the hydrostatic pressures.

Network approach for arteriolar and venular trees
The arteriolar and venular vessels are treated as a network of interconnected tubes where blood flow is described by a succession of linear relationships between the flow rate and the pressure drop [1, 14-16, 18, 20, 31]. The flow q αβ in each vessel connecting vertices α and β is written as with π α , π β the pressures at the vertices and G αβ the vessel conductance defined as with d αβ the vessel mean diameter and l αβ the vessel length. The apparent viscosity m app ab accounts for the emergence of complex microvascular flow structures (e.g. cell free layer) an their impact on average viscous dissipation at vessel scale [32]. In other words, this effective property accounts for the Fåhraeus-Lindqvist effect [33], i.e., the dependence of viscosity upon vessel diameter d αβ and discharge hematocrit H αβ . Although theoretical approaches can be used to derive expressions for the apparent viscosity (e.g. [34,35]), we use the in vivo viscosity law derived by Pries et al. in [36], where μ p represents the viscosity of the plasma, m 0:45 ab is the apparent viscosity of blood for a discharge hematocrit of 0.45 m 0:45 ab ¼ 6e À 0:085d ab þ 3:2 À 2:44e and C is a coefficient describing the dependence upon hematocrit C ¼ ð0:8 þ e À 0:075d ab Þ À 1 þ 1 1 þ 10 À 11 d 12 Further, mass balance at each inner vertex α-a vertex where no boundary condition is imposed-reads X b2N a;in where N a;in represents the set of inner neighbouring vertices that are connected to α and q s is a flow rate source term that can either describe boundary conditions for arteriolar inlets (q s > 0) and venular outlets (q s < 0) or coupling conditions for vertices connected to the homogenized capillary bed (red dots in Fig 2(a)). This latter case is detailed in Section 2.4. Alternately, when a pressure condition π s is imposed on a boundary, the corresponding vertex is called an outer vertex and Eq 6 may be rewritten X b2N a;in

Continuum approach for the capillary network
Because of its dense and space-filling structure [9], we treat the capillary bed as a porous medium where capillary vessels represent an interconnected network of pores embedded in the cerebral tissue. At the mesoscopic scale, we describe momentum transport using Darcy's law, which can be obtained from homogenization via volume averaging (see e.g. [22,23,37]) of creeping flow at microscopic scale. We have where U is the Darcy velocity (a spatially averaged velocity), P is the pressure, μ eff is the effective viscosity of blood and K eff is the permeability tensor. We assume that the capillary network is isotropic and homogeneous, so we have K eff ¼ K eff I. The effective parameters K eff and μ eff are calculated by solving microscopic-scale flow (Eqs 6 and 7) at mesoscopic scale [21]. Both coefficients are usually lumped together into a single effective resistivity K eff /μ eff of the network [14,38]. In this paper, as detailed in Appendix A, the effective permeability K eff rather represents the permeability of the network with uniform viscosity, while the effective viscosity μ eff represents the impact of the Fåhraeus-Lindqvist effect (Eq (3)). Mass balance in the porous medium is written as with q s the flow rate source term, which results from the coupling condition applied at the coupling point s located at the position x s , at any endpoint of feeding arteriolar or draining venular trees (more detail in Section 2.4). S represents the set of all sources and δ(x − x s ) the Dirac delta corresponding to the point source s. To obtain a system of equations involving only pressures, we substitute Eq 8 into Eq 9, which yields The solution of this partial differential equation is computed using a standard cell-centered finite volume (FV) approach [39,40] on a Cartesian grid. For convenience, we consider cubic cells of volume h 3 , although extension to non-cubic cells is straightforward. For cell i, the discretization reads where N i is the set containing the indices of all the neighbours to cell i, P i and P j are the approximated pressures of cells i and j, S i is the set of sources in cell i and q s is the flow rate source term. The quantified pressure P defined at the center of a cell is a local averaging of the pressure field. In the visualizations of this paper, we will consider a piecewise constant reconstruction of the pressure in the continuum, although a more precise reconstruction can be obtained by linear interpolation.

Multiscale coupling condition to link arteriolar and venular trees to the capillary bed
So far, we have described two separate frameworks on distinct scales to represent momentum transport in brain microcirculation: a network approach for arteriolar and venular trees and a homogenized continuum representation for the capillary network. Each of these descriptions is self-sufficient when boundary conditions and source terms are fixed for each separately. For instance, if q s and π s are known in the network approach, we can assemble the set of linear equations Eqs 6 and 7 into a sparse matrix plus a right-hand side and solve this linear system to obtain the pressure at each inner vertex and thus the flow rate in each vessel. However, some of the boundary conditions are in fact not known, but rather correspond to coupling elements with the continuum description. In this case, a coupling condition must be developed to express the relationship between the pressure π s at the coupled endpoint of the arteriolar or venular vessel (Eq 7) and the pressures of the local finite volume cells.
Here, we present an approach to couple these two frameworks, which is inspired from that of Peaceman [41], who developed a "well model" for petroleum engineering. In these applications, a wellbore model is coupled with a continuum reservoir description, in which variations at the scale of the underlying structure are not described. The main idea is thus to represent the strong gradients that occur in the vicinity of the wellbore, i.e., within a length scale much smaller than the size of a FV cell, by a singular pressure drop term embedded in the coupling condition. To this end, analytical expressions of the pressure field in the vicinity of the coupling points are used to derive relationships between variables of both frameworks.
In our approach, the pressure field in the neighbourhood of a coupling point is decomposed into two analytical approximations. In the closest region O s lin to the coupling point, a linear variation captures the network structure at the scale of the capillary vessels connected to this point. In O s sph , up to two FV cells away from the coupling point, a spherical solution of the continuum equations then described the variation at the scale of h, the side of a FV cell. The notations used, the relevant spatial domains and the modelling principles are schematized in Fig 3. A detailed nomenclature of the domains is also provided in Appendix D.
The main advantage of the resulting coupling condition is that it expresses a linear relationship between the pressure at the coupled endpoint of the arteriolar or venular vessel and the pressures of the local finite volume cells. Thus it can be applied to each coupling point independently. Since each coupling is blind to the others, we focus on a single coupling point in this section to explain the development of the coupling model.
Let us consider a single coupling point s, at any arteriolar or venular coupled endpoint. The local flow rate q s corresponds to the coupling term that appears in both Eq 6 for the network approach and Eq 11 for the continuum representation. The coupling condition aims at describing the pressure drop between the pressure π s at the endpoint of the arteriolar or venular vessel, and the cell pressures O s FV;neigh represents the set of cells that surround the coupled one (see Fig 3( We first present an analytical solution of Darcy's law (Eq 10) close to the coupling point, in the neighbourhood O s sph (Fig 3(b)). Assuming a spherical symmetry of the pressure field for a single coupling, we have with π ref a reference pressure defined at an arbitrary distance l ref to the source, at microscopic scale, and r ! l ref the distance to the coupling point. This solution yields a 1r-decrease of the pressure, as displayed in blue on Fig 3(d).
In order to determine π ref and l ref , a second analytical approximation of the pressure field is introduced for r l ref at microscopic scale. To this end, we extend the network approach (Eq 6) to the capillaries O cap connected to the point s (see Fig 3(c)). Replacing the subscripts β that refer to inner vertices of the arteriolar and venular trees, in Eq 6, by subscripts γ referring to inner vertices of the capillary bed, this reads X From this equation, we derive a linear pressure drop, as displayed in red on Fig 3(d), with the weighted average pressure and This expression is valid over a length scale l Γ , which defines the domain O s lin in Fig 3(c). Eq 14 can be rewritten as the pointwise formulation We then match both analytical approximations by imposing π ref = π Γ in Eq 12. This yields with The computation of l Γ is detailed in Appendix B. At this stage, the pressure field is defined at any point of R 3 in the neighbourhood of the coupling point s. Let r be the distance from a given point of the neighbourhood to the coupling point s. From Eqs 17 and 18, we have as displayed on the right-hand side of Fig 3(d).
The coupling condition is completed by assuming the existence of a domain where both the spherical approximation (Eq 20b and right part of Fig 3(d)) and the FV formulation (Eq 11 and left part of Fig 3(d)) are valid, so that the pressure of a FV cell of this domain matches with the pointwise spherical expression applied to the distance between the coupling point and the center of the cell (e.g. cell pressure P k in Fig 3(d)). This condition is based on the key assumption that h is sufficiently larger than l cap . We further consider that the matching occurs only in the most distant cells of O s sph , " O s sph . Thus, in any cell k in " O s sph , the cell pressure is expressed as where d k,s is the distance between the coupling point and the center of cell k.
Combining Eqs 14 and 21, we can obtain an expression of q s as a function of π s and ðP k Þ k2 " where P j k represents the pressure of the unique neighbouring cell j k of k, which belongs to " O s FV;neigh (see Fig 3(b)). By combining Eqs 14, 21 and 22, we finally obtain the coupling expression where f is the following linear function with Card(Á) representing the cardinal number of a given set and Recall that the FV cells are assumed to be cubic here, although the extension to parallelepipedic cells is straightforward. The development above is valid for both centered and off-centered couplings, as in Fig 3(e). However, the FV scheme (Eq 11) is not sensitive to the position of the coupling point in a single source cell. Thus, the resulting FV pressure field is independent of the source position in a coupled cell. To properly take into account the off-centering of a coupling point s in the FV part of the hybrid approach, the coupling term q s is redistributed on each cell i of the neighbourhood O s FV;neigh with new partial source terms defined by according to a partition coefficient τ i,s where v i is the volume of the intersection between the cell i and a fictitious mesh cell centered on the coupling point (see S1 Fig in Supporting Information). For a given source s, this definition implies that ∑ i τ i,s = 1, yielding ∑ i q i,s = q s , and that, among the 27 cells of O s FV;neigh no more than 8 verify τ i,s 6 ¼ 0, and q i,s 6 ¼ 0. From the point of view of the FV scheme (Eq 11), a given cell i of the FV domain is impacted by any coupling located in its neighbourhood " N i , which includes the cells sharing a single corner with i. Thus, Eq 11 is can be rewritten as where Sð " N i Þ is the set of coupling points in " N i . Combining Eqs 24 and 26-28, Eq 28 finally reads

Matrix structure and numerical implementation
Due to the linearity of the hybrid model, we can compute the flow in the whole domain by solving a single linear system, assembled from the equations describing the FV (Eq 28), the network (Eqs 6 and 7) and the coupling (Eqs 23-25). This yields a sparse matrix which is composed of seven parts: two symmetric sparse blocks for the description of the FV and the network, of size (number of cells) 2 and (number of inner vertices) 2 , respectively; a diagonal block, of size (number of couplings) 2 ; and four supplementary off-diagonal blocks for the description of the relationships between the couplings points and both the network and the continuum (see more detail on the matrix structure in Appendix C). The right-hand side describes the boundary conditions imposed at the inlets and outlets of the arteriolar and venular trees, and at the six faces of the FV domain. The sparsity of the linear system allows us to take advantage of the efficient and scalable numerical tools in the PETSc library [42]. We use a block-Jacobi preconditioner, which is very well-suited to parallel computations, and the system is solved via the bi-conjugate gradient iterative method. The code was developed from scratch in C++ and is fully parallelized using MPI, PETSc [43] and the ParMETIS library for parallel graph partitioning of the network [44]. The code is designed for High Performance Computing (HPC) and runs on EOS supercomputer (CALMIP, France). Regular automatic testing, using CxxTest (http://cxxtest.com/), ensures reliability of the code [45].
Alternately, this code can be used to run computations where all the vessels, including the capillaries, are modelled by the network approach. These complete network (CN) simulations will be used in Section 3 to generate reference fields enabling the evaluation of the accuracy of the hybrid approach.

Architecture of capillary networks and modelling configurations
With a view to model validation, calculations are performed here for two types of model capillary networks, whose names were chosen according to graph theory [46]: (1) 6-regular networks, i.e., where each vertex is connected to 6 neighbours; and (2) 3-regular networks, which comply to the physiological connectivity in healthy capillary networks (see Fig 4). The use of synthetic capillary networks enables both to solve the problem of anatomical data availability at this scale and to verify all the modelling assumptions made in the present paper (isotropy and uniform permeability). Furthermore, using anatomical capillary datasets would generate uncertainties in the computation of the effective parameters that characterize the capillary bed in the continuum representation. This would add errors in the hybrid resolution that would be indiscernible from errors induced by the continuous/discrete coupling model, going against our wish to investigate the intrinsic impact of the multiscale coupling model on large-scale flow fields.
In both 6-and 3-regular networks, the capillaries have uniform lengths, l cap = 50 μm in 6-regular networks and l cap = 30.62 μm in 3-regular networks, so that the two network types have similar vascular densities. Capillary diameters are either uniform (d cap = 8 μm) or distributed following a Gaussian distribution with a mean diameter of 5.91 μm and a standard deviation of 1.30 μm, which is characteristic of human capillary networks [8]. The corresponding effective parameters K eff and μ eff for the continuum representation are presented in Table 1 (more detail in Appendix A).
We will first consider basic configurations that involve less than two couplings. For a single coupling, an arteriole penetrates into a cubic 6-regular capillary bed of size (2.25 mm) 3 containing about 284000 vessels. For the case of two coupling endpoints, an arteriole and a venule plunge into a parallelepipedic capillary bed of size 5.5 mm × 2.25 mm × 2.25 mm containing  Table 1. Effective parameters for the continuum representation of the model capillary networks. This table summarizes the values of effective permeability K eff and viscosity μ eff used in this paper (Section 2.6). Details about the computation of these coefficients are provided in Appendix A. 557000 vessels. Because we want to investigate the impact of the cell size h on the accuracy of the model, we choose to ensure mesh convergence by considering large capillary domains. Thus we also avoid the potential interference of boundaries with the reconstruction of the pressure and flow rate fields. The arteriolar and venular vessels are l AV = 200 μm long and d AV = 15 μm in diameter. The boundary conditions are illustrated in Fig 5. For consistency, the boundary conditions imposed at the boundaries of the capillary bed in the CN simulations are similar to those imposed at the faces of the continuum in the hybrid approach. In the case with two couplings (Fig 5(b)), we impose a constant pressure drop δP = 7000 Pa between the arteriolar inlet and the venular outlet [15,47]; impermeability (zero flow) on the top and bottom of the capillary network; and periodic conditions in the two other directions. The condition of impermeability imposed at the top of the domain is physiologically relevant to describe the surface of the cortex. The same condition is usually imposed at the bottom, e.g. [14,15]. For a single coupling (Fig 5(a)), boundary conditions are similar with a zero pressure condition at the bottom of the capillary bed, inducing a global pressure drop δP = 10000 Pa.
Second, we will consider a more realistic case, which involves the anatomical arteriolar and venular trees presented in Fig 2(a). This dataset has been previously obtained by Cassot et al. [8] from sections of a human brain, from the Duvernoy collection [7]. The dataset supporting this article is provided in a supplementary data file. It is composed of a venular tree in between two arteriolar trees, containing a total of 563 vessels and displaying 274 coupling points with a realistic spatial distribution, as displayed on Fig 2(b). Arteriolar and venular vessel diameters are in the range d AV 2 [7.5 μm; 42.6 μm]. They are connected to a 3-regular capillary network of size 4 mm × 3 mm × 1.5 mm containing 1.6 million vessels. The boundary conditions are identical to the ones defined above for two couplings, as shown in Fig 5(b), except that the pressure drop of 7 000 Pa is imposed between the two inlets of the arteriolar trees and the outlet of the venular tree.

Error metrics
To evaluate the accuracy of the model, we compare results of the hybrid approach with CN simulations, where blood flow in all vessels is described by the network approach (Eqs 6 and 7). We consider the four domains of interest illustrated in Fig 3( where p represents the pressure value, the subscripts H and CN refer to the hybrid and CN approaches, and the subscript ref refers to the reference value used for normalization. For all component ω of O q , the local flow rate error reads where q represents the flow rate value. The global errors characterizing O are further defined as and For the study of errors distribution, we use the local errors in Eqs 30 and 31. To avoid underestimation, these errors are normalized by local reference values, so that

Validation of the hybrid approach for basic configurations
Here, our goal is to assess the robustness of the coupling model and the relevance of strategies that have been adopted during its development. The first test case analyzes the impact of the size of FV cells on the accuracy of the model (Fig 6i(a)). Three other configurations test different kinds of perturbations that typically occur in brain microcirculation: off-centering (Fig 6ii  (a)), effects of boundary conditions (Fig 6iii(a)) and coupling points coming close together (Fig 6iv(a)). The accuracy of the model is determined by comparing the hybrid and CN approaches, using the global error metrics, Eqs 32 and 33, for the four domains of interest illustrated in Fig 3(e). While these errors may not be as representative as local errors, this choice enables us to compare the different graphs displayed in Fig 6 and to conclude about the importance of the impact of each perturbation.
An important hypothesis of the coupling model is that the size h of a mesh cell is sufficiently larger than the characteristic length l cap of the capillaries-so that Eq 21 is valid. The most sensitive area to this condition is a priori the distant neighbourhood O s;ε sph (Fig 3(e)), which corresponds to the intermediate zone between the analytical approximations and the FV formulation. In the first test case (Fig 6i), we focus on identifying the threshold value of h/l cap for which the FV cells are large enough. As expected, the errors ε p  (Fig 6i(b) and 6i(d)). For comparison, the pressure and flow rate errors have also been computed when imposing a simple condition of pressure continuity at the interface between the arteriolar vessel and the continuum (see S2 Fig in Supporting Information). In this case, all errors significantly increase with the ratio h/l cap . In particular, the pressure and flow rate errors at the coupling point are * 29% and * 150%, respectively, for h/l cap = 4. They continue to increase for larger values of h/l cap , reaching * 45% and * 250%, respectively, for h/l cap = 9. This behaviour highlights the difference in the physical meanings between the pressure of the coupled FV cell and the pressure at the endpoint of the arteriolar vessel, which is critical for large values of h/l cap . By contrast, the multiscale coupling model included in our approach takes into account the upscaling transition between the arteriolar vessel and the continuum. Moreover, other geometrical perturbations such as nonuniform capillary diameters distribution, a non-vascular zone around the arteriolar vessel representing the Virshow-Robin space [7] or a different architecture of the capillary bed have no influence on the threshold value of h/l cap (S3 Fig in Supporting Information). In the remainder of this paper, we therefore always use h/l cap ! 4.
As displayed in Fig 6i(c)), for a ratio h/l cap = 5, i.e. h = 250 μm, the pressure errors in the capillary medium O ε O FV are lower than 1% everywhere when using the multiscale coupling model, showing its accuracy. Furthermore, the reconstruction of the pressure field using the analytical approximations (Eqs 14 and 12) in the vicinity O s FV;neigh of the coupling point is in very good agreement with the reference CN pressure (Fig 7i). As expected, the FV representation is accurate beyond two cells away from the coupling but does not describe the strong gradients in O s FV;neigh , underestimating the pressure in the coupled cell. The coupling model correctly captures these gradients since the pressures π CN and π H at the coupling point are in very good agreement, with a global error of 0.7%. By comparison, the difference between these two values is 100 times larger when using a simple condition of pressure continuity, with a global error of 33% (S2 and S4 Figs in Supporting Information). Moreover, because the additional pressure drop building up around the coupling point is not accounted for, the pressure at the arteriolar endpoint is highly underestimated (π/π ref = 0.39 vs. 1.05 in the CN approach, see S4 Fig). Thus, this local singular pressure drop significantly contributes to the global resistance of the capillary bed, which has been recently shown to dominate the overall intracortical resistance in both human and mice [31,48]. In other words, this microscopic scale effect, intrinsically linked to the transition between the tree-like arterioles and the mesh-like capillaries, has a huge impact on the pressure field at macroscopic scale.
We now ask whether the off-centering of the arteriolar or venular endpoints in the finitevolume cell affects the accuracy of the result. Off-centering of coupling points can indeed induce asymmetries of the pressure field with regard to the FV grid. To deal with this, we introduce partial sources (Eq 28) to distribute the source term among the cells affected by the well model and compare the results with computations where the flow rate is not distributed. In Fig 6ii, we analyze the impact on the pressure and flow rate errors (Fig 6ii(b) and 6ii(d)) by progressively increasing the distance d off to the FV cell center in the range d off = l cap 2 0:09; ffi ffi  6ii). By reconstructing the analytical approximation of the pressure field in the vicinity of the coupling (Fig 7ii), we see that the coupling model without distribution of q s cannot capture the asymmetry of the pressure field. This induces slight over-and underestimations on each side of the coupling point that persist even far away from the source. On the contrary, the reconstructions with q s distributed are in very good agreement with the reference CN computations.
Next, we assess the influence of boundary conditions in Fig 6iii. The coupling point is initially positioned at the center of the domain (d boundary /h = 4.5) and progressively brought close to a boundary (d boundary /h = 0.5) where impermeability conditions have been imposed. Similarly to the first case, h = 250 μm here for the purpose of working with centered couplings and avoid potential errors induced by off-centering. Both pressure and flow rate errors (Fig 6iii(b)  and 6iii(d)) show that our model is robust provided that at least one cell separates the coupling point from the domain boundary, ensuring that the neighbourhood O s FV;neigh of the coupling point does not intersect it. While the model is less accurate when the coupling occurs in the boundary cell, the increase in errors is small (up to 2.5% in pressure and 6% in flow rate).
Finally, we study in Fig 6iv the impact of the distance between an arteriolar coupling and a venular coupling on results, as the density of coupling points can be locally very large in the human cortex (Fig 2). Initially separated by 9 cells of side h = 250 μm (d AV /h = 9), the two couplings are brought closer together until they are in neighbouring cells (d AV /h = 1). In Fig 6iv  (b) and 6iv(d), both pressure and flow rate errors show that the hybrid approach is in very good agreement with the CN reference solution, even for short distances that break the spherical symmetry of the pressure field. For the ratio d AV /h = 1, the strong decrease of errors can be explained by the fact that each coupled cell belongs to the neighbourhood O s FV;neigh of the other one. In this case, a shunt appears at the scale of the FV cells and the pressure field is uniform outside the coupling region.
In summary, the hybrid approach requires a coupling model much more complex than a simple condition of pressure continuity to capture the strong gradients around coupling points. The coupling condition presented in this paper fills this gap since it is robust to the geometrical perturbations that typically occur in the brain microcirculation. For more than two coupling points, we also verified that the model is accurate in similar idealized configurations for up to 10 couplings occurring in the same cell ( S5 Fig in Supporting Information). For comparison, the results are summarized in Fig 8. The values reported correspond to the maximal errors computed for h/l cap > 4 (orange dashed lines displayed on each previous graphs- Fig 6  and S5 Fig). With regard to pressure errors, the hybrid approach always shows very good agreement with the CN approach, with no error larger than 3%. The largest errors are observed for the configuration involving multiple couplings in one FV cell, suggesting that the accuracy of the reconstruction of the pressure field in the continuum is limited by the density of couplings per cell. This trend is confirmed when considering flow rate errors in Fig 8. In the next section, we explore this in more detail through a more realistic case representative of the architecture of the human vascular network, where up to 5 couplings occur in the same cell.

Application to larger realistic networks involving anatomical arteriolar and venular trees and synthetic capillary networks
The errors presented in Fig 9 refers to a more realistic configuration involving the three anatomically accurate arteriolar and venular trees displayed in Fig 2(a), which plunge into a 3-regular capillary network with Gaussian diameters distribution and uniform length l cap = 30.62 μm. The corresponding continuum is discretized in FV cells of side h = 122.5 μm, i.e. h/l cap = 4, ensuring the validity of the coupling model. In this case, the global flow rate error, Eq 33, computed at coupling points is very low (<0.26%). Consequently, this also ensures a very good accuracy in the continuum. To go further, we analyze the distribution of local errors, Eqs 30 and 31, in the arteriolar and venular trees.
The pressure error field (Fig 9(a)) shows that the hybrid approach remains accurate, with no error greater than 6%. Flow rate is much more sensitive, with errors ranging from 0.001% to 100%. These large values can be explained by the fact that pressure and flow rate fields are highly correlated (Eq 1) in the sense that slight pressure perturbations at the endpoints of a given vessel can induce very high flow rate error, especially for low pressure drops (low flow rates). In spite of this, considering both pressure and flow rate errors in Fig 9(b), we clearly observe the benefit of using the coupling model. This is confirmed when considering absolute values instead of relative errors (S6 Fig in Supporting Information). Moreover, the large errors highlighted in vessels with low flow rates have a minor impact on the spatial distribution of the source flow rates. This can be verified by comparing the cumulative flow rates, e.g. in the left arteriole (Fig 10(a)), deduced from the hybrid approach and CN computations. For that purpose, the source flow rates q s are first sorted in ascending order and the set of cumulated flow rates is defined for n in {1, . . ., Card(S)}, where S is the set of sources, by This means that Q CardðSÞ s must be equal to the flow rate q A,l injected at the entry of the left arteriole. Finally, the results are normalized by q A,l . They show very good agreement with CN computations, which is not the case when considering a simple condition of pressure continuity in the hybrid approach (Fig 10(a)). Moreover, probability density functions of the flow rate in the network vessels (Fig 10(b)) also show very good accuracy as the Gaussian fits overlap, contrary to those obtained without coupling. These results are particularly interesting because the standard complete network approach has been validated only in a statistical sense, rather than a vessel-to-vessel basis, by comparison with in vivo data [20].
Finally, the results presented in this section confirm that the coupling model is very accurate, even for multiple couplings.

Computational performance
In order to estimate the performance gain of the parallel design of the code, blood flow was computed in a network containing more than 10 millions of vertices. This lattice was obtained by periodically duplicating the pattern of three arteriolar and venular trees (Fig 9). The final network contains 10.1 million of vertices in the CN version: 50000 arteriolar and venular vertices and 10.05 million capillary vertices, these latter being homogenized into 51000 mesh cells in the hybrid approach. The computation was run on 20 cores Intel(R) Xeon(R) CPU E5-2680 v2 @ 2.80GHz of the EOS supercomputer (CALMIP, France). While the CN approach takes about 6 min, the hybrid approach takes only 1 s, implying a ×360 acceleration. This important gain in computation time is explained by two facts: (1) the number of unknowns is reduced by a factor of *200 in the hybrid approach; (2) the unstructured architecture of the capillary network is replaced by a structured grid, which requires much less numerical treatment. By using a simple linear projection of this result, the flow computation in the whole human cortex (*1000 cm 3 [49] with *10 000 vessels per mm 3 [8], i.e. 10 billions of vessels) would take about four days with the CN approach while it would take *15 minutes with the hybrid approach. This demonstrates that the combination of a HPC code with the hybrid approach is a powerful tool that pushes the current limits of computational performance for simulating blood flow in the brain microcirculation.

Comparison with previous work
Now that we have established the relevance and the numerical strength of our approach, in this section we emphasize the main methodological differences between our model and previous work.
A few authors have introduced hybrid approaches where the arteriolar and venular trees are represented by a discrete model and the capillary network as a continuum [27]. However, the need for a proper coupling condition at the interface between the two frameworks has never been discussed, even in the context of cardiac microcirculation where multiscale computational approaches have been used for a longer time [25,26,38]. Reichold et al. [14] got around this difficulty by replacing the capillary mesh by a coarser network with identical resistivity. The advantage of such method is to decrease the number of degrees of freedom of the problem, still using a complete network approach. However, this coarsening results in a larger capillary length. Thus the 1r pressure relaxation highlighted in this paper around the couplings is likely shifted to larger length scales and the strong pressure gradients around the couplings at microscopic scale may not be properly captured.
One simple alternative idea to capture these strong gradients would be to sufficiently refine the finite volume mesh close to the coupling points. However, the solution of the diffusion equation in the FV scheme (Eq 8) would impose a 1r-decrease of the pressure in the vicinity of the coupling point, overlooking the linear variation at microscopic scale, below l cap . Furthermore, even if this were physically accurate, couplings are relatively dense (Fig 2(a)), implying that such a refinement would be computationally limiting. This might still be tractable on a larger number of cores, but only to calculate blood flow in moderately large portions of the brain, without the possibility to further consider mass transport modelling or nonlinear effects in the blood rheology.
In the present hybrid approach, we used analytical solutions close to the sources to derive a proper coupling model. This idea was introduced by Peaceman in the context of petroleum engineering [41], in order to couple wellbores with two-dimensional finite volume models representing the flow in an oil reservoir. It was later extended to deal with more complex situations, e.g. anisotropic continuum, non-Darcy flow, three-dimensional models, finite element methods [50][51][52]. However, because of differences in scales between oil reservoirs and the human cortex (see S7 Fig in Supporting Information), important modifications to the original well models were made. First, the sources and sinks terms are lineic in Peaceman's approach, while our conceptual model involves point terms. Thus analytical solutions in classical well models exhibit cylindrical symmetry, while the symmetry is spherical here. More importantly, the reference radius and wellbore pressure used for the analytical solution in Peaceman's model [41] are straightforward because the well diameter is much larger than the pores [53]. In microvascular networks, the diameter of a coupled arteriole or venule is * 5 times smaller than the characteristic length of a capillary vessel [9]. Consequently, our model required an additional analytical approximation of the pressure drop in the direct vicinity of couplings, Eq 14, to capture microscopic effects. In addition, the FV cells in an oil reservoir are a thousand times larger than the well diameter [53], while they are only ten times larger than the characteristic length of a capillary in the present work. We therefore expanded the coupling condition to a larger neighbourhood O s FV;neigh around coupling points. In doing so, the model yields an additional local linear system per coupling, which involves the coupled cell and its 26 neighbours. Even if this strategy is constraining, it is a reliable way of improving the accuracy of the model [54], which is similar to extending the stencil in computational methods. Moreover, it enables us to deal with off-centered couplings by distributing the coupling condition among these particular cells, Eqs 27 and 28 (see similar strategies in [55]).

Conclusions and perspectives
In this paper, we presented a hybrid multiscale approach to simulate blood flow in large volumes of the brain microcirculation. This approach is hybrid in the sense that it couples a network representation of the arteriolar and venular trees to a homogenized continuum representing the capillary bed. We showed that the coupling of the two frameworks is nontrivial since the effects induced by this coupling at microscopic scale have a huge impact on the flow at macroscopic scale. As a consequence, a simple condition of pressure continuity yields large errors in both pressure and flow rate fields. To deal with this issue, we introduced a coupling model that relies on analytical approximations of the pressure field in the vicinity of coupling points. Such a theoretical formulation takes into account the details of the capillary network architecture in the close neighbourhood of the arteriolar and venular endpoints (number, length and diameter of capillaries connected to each coupling point).
By comparison with a complete network approach, we showed that the resulting hybrid approach is accurate in describing the pressure field in various simple test configurations, which involve model capillary lattices with strong topological and morphological differences (6-regular vs. 3-regular networks, networks with uniform or distributed diameters, impact of Virshow-Robin space). We also demonstrated the validity of the hybrid approach when interactions take place between several coupling points (either in the same grid cell or not), or with the domain boundaries. We finally validated this approach on a more complex configuration, combining a few arteriolar and venular trees derived from human anatomical data, with a model capillary network. This last configuration displays 274 coupling points, combining all the geometrical perturbations that typically occur in brain microcirculation (e.g. off-centering, close couplings). Finally, we were able to conclude that the accuracy of the hybrid approach in the large tree-like vessels is similar, on average, to that of a complete network approach.
To go one step further and use anatomical data for the representation of the capillary bed in the future, the model might require several extensions. For instance, the stratification of the vascular density in the depth of the cortex [56] and possible anisotropies [7] may induce strong heterogeneities of capillary networks. To describe this accurately, analytical approximations of the pressure field that capture non-spherical symmetries in the vicinity of couplings may be necessary (see, e.g. [57] for similar problems in petroleum engineering). Moreover, in a perspective of extension to the simulation of mass transfers, it is noteworthy that Darcy's law is similar to a diffusion equation in pressure. Introducing a hybrid approach for mass transfers will thus involve issues similar to the ones described in the present paper. In the particular case of oxygen transfers, hematocrit heterogeneities induced by the non-proportional distribution of the red blood cells at microvascular bifurcations, i.e. phase separation effect, should also be considered. While this problem has already been investigated in microvascular networks (e.g. [15,18,58]), the introduction of such local bi-phasic phenomena in the continuous description of the capillary bed is still a fully open question. Since hematocrit heterogeneities induce non-linear effects, this will likely imply solving the linear system introduced in this paper in an iterative fashion. Thus, the improvement in computational performance will be even more critical for whole-brain computations, not to mention the additional resources needed to solve for the mass transfer problem.
In this perspective, the present hybrid approach has been developed for HPC. With regard to blood flow, this will enable (1) the validation of the hybrid approach by comparison with CN computations for a cortical volume of up to 1 cm 3 (about ten million of vessels), reaching the computational limit of the CN approach; (2) the simulation of blood flow in much larger volumes-with the ultimate goal of modelling the whole human cortex-via the hybrid approach, assuming that sufficient anatomical data becomes available for arteriolar and venular networks [59]. Such simulations would represent a valuable tool for studying flow variations induced by vascular modifications-vascular occlusions or rarefaction of the capillary vessels observed in Alzheimer's disease.

Appendix A Effective properties of the capillary bed
In this appendix, we present how the effective permeability K eff and viscosity μ eff are calculated in 6-and 3-regular capillary networks displayed in Fig 4. Both coefficients are usually lumped together by considering a single effective resistivity K eff /μ eff of the network. Here, the effective permeability K eff represents the permeability of the network with uniform viscosity, while the effective viscosity μ eff represents the impact of the Fåhraeus-Lindqvist effect (Eq 3).
To compute these effective parameters, the flow is first solved at microscopic-scale in cubic capillary networks of side L, using the network approach (Eqs 6 and 7). The results are then averaged at mesoscopic scale and identified to Darcy's law where subscript L refers to the network side, Δ L P is the imposed pressure drop and Q L is the resulting global flow rate. The side L of the domain is increased until the effective coefficients do not depend on it anymore, meaning that we have reached the mesoscopic scale at which the capillary bed can be considered as a continuum. The side L that corresponds to the convergence of the effective coefficients defines what is commonly called a Representative Elementary Volume (REV) [21,23]. The characteristics of the capillary networks treated and the final results are summarized in Table 1. In the 6-regular network, we can theoretically derive the corresponding value of K eff . In a cubic domain of side L, this network is indeed similar to N × N parallel vessels of length L, where N represents the number of vessels normal to the direction of interest. We have with the conductance corresponding to a vessel of length L Noting that L/N = l cap here, we deduce from Eqs 35 and 36 a value of K eff that does not depend on L This is valid for L > l cap . Here, K eff ' 4.02 × 10 −14 m 2 in each direction, with a REV side equal to l cap .
No theoretical formula can be derived to calculate the corresponding permeability K eff in the 3-regular network. The REV side and the associated effective parameters are determined by increasing the side of the network L until a convergence is observed. As displayed in Fig 11, K eff L finally reaches a plateau corresponding to the converged effective permeability K eff . The beginning of the plateau corresponds to the side of the REV. We set the convergence criteria at |K(L n ) − K(L n−1 )|/K(L n ) < 0.01, where L n−1 represents the domain side just below a given sideL n . Here, the REV side is about 300 μm and the corresponding permeability value is K eff ' 6.61 × 10 −14 m 2 . Appendix A.2 Distributed capillary diameters. If the diameter is non-uniform, it follows the Gaussian law detailed in Section 2.6. Since the capillary network is statistically generatedso the diameters distribution follows the Gaussian law, we study the statistical properties of the resulting permeability and viscosity, following the three steps below (results are presented in Fig 12 for both 6-and 3-regular networks): 1. In order to quantify the impact of a diameter distribution on the effective parameters K eff and μ eff , we first compute the reference values K ref and μ ref in the case of uniform capillary diameters d const = 5.91 μm (Fig 12(a)), which corresponds to the mean diameter of the Gaussian distribution. From the in vivo viscosity law, we have μ ref ' 5.95 × 10 -3 Pa Á s for both networks. The reference permeability K ref is then deduced from Eq 35 (Fig 12(a)).  Fig 12(b). We establish the convergence for a standard deviation under 0.1. For the 6-regular   Fig 12(c). For the 6-regular network, μ eff ' 5.88 × 10 -3 Pa Á s at the REV side of 200 μm. For 3-regular network, μ eff ' 6.15 × 10 -3 Pa Á s at the REV side of 300 μm.
As expected, we also deduce that the two networks are isotropic because the results overlap in the three directions of space (Figs 11 and 12). Indeed, all networks used in the present paper are periodic in the x, y and z directions, with symmetry planes along the three axis. Thus, according to [21], even if the pressure gradient has not the same principal orientation as the segments of the network, Darcy's equation still holds. In other words, while segments have preferential orientations at microscopic scale, the medium is isotropic at mesoscopic scale.
Multiplying by G sγ and summing over V s , we obtain which also can be written as Finally, by identifying π Γ and l Γ , this equation corresponds to Eq 21 with Appendix C Matrix structure The matrix describing the linear system for the hybrid approach is a sparse matrix, which is composed of seven parts displayed on Using the hybrid approach, we can compute blood flow by solving a single linear system, assembled from the equations describing the finite volume discretization in the continuum (FV, in green), the network approach in the arteriolar and venular trees (AV, in black) and the coupling model at the interface of the two frameworks (coupling, in red). This yields a sparse matrix that is schematized here. Plain lines and dots represent non-zero values. The components (P i ) i , (π α ) α and (π s ) s are the unknowns of the linear system, corresponding to FV, AV and coupling parts, respectively. https://doi.org/10.1371/journal.pone.0189474.g013 Part c b : the coupling part is also completed by an off-diagonal block of size (number of cou-plings×number of inner vertices). The components correspond to the connection of each coupled vertex to its unique neighbouring vertex in the network part.
The right-hand side, which describes the boundary conditions imposed at the inlets and outlets of the arteriolar and venular trees, is finally very sparse.

Appendix D Nomenclature
In this section, a detailed nomenclature is provided about the description of the domains used in the main paper. The domains are sorted according to the elements that they contain: network domains contain vessels, finite volume domains contain discretization cells, and continuous domains are subdomains of R 3 . Each of the domain detailed below is assigned to one of the following classes: main structural component, domain involved in the development of the coupling model, validity domains of the analytical approximations, and domain dedicated to the computation of errors.

Appendix D.1 Network domains.
• O AV is the set of network vessels in the arteriolar and venular trees (main structural component); • s is a coupling point, at the interface between O AV and O FV (main structural component); • O cap is the set of capillary vessels connected to a coupling point (main structural component).

Appendix D.2 Finite volume domains.
• O FV is the set of cells used to discretize the continuum (main structural component);  In complete network version, arteriolar and venular trees are connected to 6-or 3-regular capillary networks. In the hybrid approach, they plunge into an equivalent porous medium which is discretized in finite volume cells of size h. (a) In configuration i, a single coupling is imposed at the center of the coupled cell. We then vary h/l cap ratio where l cap = 50 μm represents the characteristic length of a capillary vessel. In the following configurations, the same study has been led by taking into account: ii. the Virshow-Robin space (nonvascularization around arteriole [7]); iii. the diameters distribution; and iv. the physiological 3-regularity of the capillary network.  Fig 6i(a) in the main paper. The reference pressure field from CN approach is represented by orange crosses and green steps represent cell pressures of the FV representation. The pressures π H and π CN at the endpoint of the coupled vessel, for both approaches, are indicated by tick red and orange dashes, respectively. Here π ref = 5000Pa is a reference pressure value used for nondimensionalization. (EPS) S5 Fig. Global pressure and flow rate errors for multiple couplings occuring in the same FV cell. From 2 to 10 arteriolar and venular vessels are randomly coupled to a 6-regular capillary network, in such a way that they are located in the same FV cell in the FV domain of the equivalent hybrid configuration. According to physiological statistics (Fig 2(b)), a minimal distance of 100 μm is imposed between the arteriolar inlets and venular outlets.