Maxime Berg

Yohan Davit

Michel Quintard

Sylvie Lorthois

Modelling solute transport in the brain microcirculation: is it really well mixed inside the blood vessels?

Keywords:

Most network models describing solute transport in the brain microcirculation use the well-mixed hypothesis and assume that radial gradients inside the blood vessels are negligible. Recent experimental data suggest that these gradients, which may result from heterogeneities in the velocity field or consumption in the tissue, may in fact be important. Here, we study the validity of the well-mixed hypothesis in network models of solute transport using theoretical and computational approaches. We focus on regimes of weak coupling where the transport problem inside the vasculature is independent of the concentration field in the tissue. In these regimes, the boundary condition between vessels and tissue can be modelled by a Robin boundary condition. For this boundary condition and for a single cylindrical capillary, we derive a one-dimensional cross-section average transport problem with dispersion and exchange coefficients capturing the effects of radial gradients. We then extend this model to a network of connected tubes and solve the problem in a complex anatomical network. By comparing with results based on the well-mixed hypothesis, we find that dispersive effects are a fundamental component of transport in transient situations with relatively rapid injections, i.e. frequencies above one Hertz. For slowly varying signals and steady states, radial gradients also significantly impact the spatial distribution of vessel/tissue exchange for molecules that easily cross the blood brain barrier. This suggests that radial gradients cannot be systematically neglected and that there is a crucial need to determine the impact of spatio-temporal heterogeneities on transport in the brain microcirculation.

Introduction

The human cerebral circulation can be decomposed into three main components: the large surface vessels (pial arteries and veins, figure 1a) that connect the brain to the heart; smaller vessels (arterioles and venules) that branch off from these large vessels and dive into the cortex (figure 1b); and the capillary bed forming a complex interconnected network that perfuses throughout the grey matter (figure 1c). Together, ....... ....... the arterioles, venules and capillary bed represent the brain microcirculation. They play a central role in cerebral homeostasis as they control exchanges between the vasculature and the brain tissue (fi gure ld). This includes the delivery of vital molecules (oxygen and nutrients) to neurons and the removal of toxic wastes (e.g. CO 2 , amyloid) through the blood-brain barrier, i.e. the semi-permeable layer of endothelial cells forming the walls of cerebral microvessels [START_REF] Abbott | Structure and function of the blood-brain barrier[END_REF]. The microcirculation also plays a key role in various pathologies, e.g. vessels that get stalled or clogged are involved in stroke and neurodegenerative diseases such as Alzheimer's disease [START_REF] Gorelick | Vascular contributions to cognitive impairment and dementia[END_REF][START_REF] Zlokovic | Neurovascular pathways to neurodegeneration in Alzheimer's disease and other disorders[END_REF][START_REF] Brundel | Cerebral microinfarcts: a systematic review of neuropathological studies[END_REF]Cruz-Hernandez et al. 2019). Accurate models of blood fl ow and molecular transport in large microvascular networks are needed in order to better understand the mechanisms associated with these pathologies, to translate results from animal models to humans and to assess the efficacy of new treatment strategies (Cruz-Hernandez et al. 2019). Such models could also improve diagnosis by enriching the intra voxel models of mass transport (see e.g. [START_REF] Barnes | Modeling the effect of intra-voxel diffusion of contrast agent on the quantitative analysis of dynamic contrast enhanced magnetic resonance imaging[END_REF] used in the interpretation of clinical imaging techniques such as perfusion or functional magnetic resonance imaging [START_REF] Holdsworth | Magnetic resonance imaging techniques: fMRI, DWI, and PWI[END_REF].

Developing such microvascular transport models is challenging. The main difficulty in doing so is to accurately capture the multiscale and multiphysics nature of the problem, which couples nonlinear effects in the rheology of blood [START_REF] Pries | Red cell distribution at microvascular bifurcations[END_REF][START_REF] Pries | Microvascular blood viscosity in vivo and the endothelial surface layer[END_REF], molecular transfers across complex biological tissues [START_REF] Abbott | Structure and function of the blood-brain barrier[END_REF][START_REF] Kutuzov | Contributions of the glycocalyx, endothelium, and extravascular compartment to the blood-brain barrier[END_REF] and effects linked to the FIGURE 2. Flow rates and Péclet numbers within a mouse microvascular network, as computed in [START_REF] Bracko | Neutrophil adhesion in brain capillaries contributes to cortical blood flow decreases and impaired memory function in a mouse model of Alzheimer's disease[END_REF]. (a) Flow rates in a network that contains more than 15 000 vessels for a volume of approximately 1 cubic millimetre. It was obtained post mortem and imaged by two photon microscopy [START_REF] Tsai | Correlations of neuronal and microvascular densities in murine cortex revealed by direct counting and colocalization of nuclei and vessels[END_REF][START_REF] Blinder | The cortical angiome: an interconnected vascular network with noncolumnar patterns of blood flow[END_REF]). The topology is described by an adjacency matrix, storing the list of all edges (e.g. 0, 1 and 2) connected to a given vertex (e.g. B), as illustrated for a single bifurcation. (b) Distribution of the radial Péclet numbers, U R/D V , based on the blood flow distribution shown in (a) for a diffusion coefficient D mol = 10 -9 m 2 s -1 . Each dot represents a different vessel and the coloured surface shows the local density of points.

geometry/topology of the vasculature [START_REF] Lorthois | Fractal analysis of vascular networks: insights from morphogenesis[END_REF]. The most widespread approach to modelling such systems, is to consider each microvessel as a tube and the microvasculature as a network of connected tubes (illustrated in figure 2(a), see also [START_REF] Pries | Blood flow in microvascular networks. Experiments and simulation[END_REF], [START_REF] Goldman | Computational modeling of oxygen transport from complex capillary networks. Relation to the microcirculation physiome[END_REF], [START_REF] Secomb | Green's function methods for analysis of oxygen delivery to tissue by microvascular networks[END_REF], [START_REF] Fang | Oxygen advection and diffusion in a three-dimensional vascular anatomical network[END_REF], [START_REF] Obrist | Red blood cell distribution in simplified capillary networks[END_REF], [START_REF] Linninger | Cerebral microcirculation and oxygen tension in the human secondary cortex[END_REF], [START_REF] Kojic | A composite smeared finite element for mass transport in capillary systems and biological tissue[END_REF], [START_REF] Peyrounette | Multiscale modelling of blood flow in cerebral microcirculation: details at capillary scale control accuracy at the level of the cortex[END_REF], and [START_REF] Sweeney | Insights into cerebral haemodynamics and oxygenation utilising in vivo mural cell imaging and mathematical modelling[END_REF]). The advantage of this representation is that it simplifies the treatment of each vessel, making it possible to model large volumes and to capture the global topology of the network. The disadvantage, however, is that its accuracy heavily relies on the quality of the models used to represent each individual tube and the couplings between tubes. This framework is well established for blood flow [START_REF] Pries | Blood flow in microvascular networks. Experiments and simulation[END_REF][START_REF] Lorthois | Simulation study of brain blood flow regulation by intra-cortical arterioles in an anatomically accurate large human vascular network. Part I. Methodology and baseline flow[END_REF][START_REF] Fry | Estimation of blood flow rates in large microvascular networks[END_REF]) and accurately describes both in vivo and in vitro data [START_REF] Pries | Microvascular blood viscosity in vivo and the endothelial surface layer[END_REF][START_REF] Bracko | Neutrophil adhesion in brain capillaries contributes to cortical blood flow decreases and impaired memory function in a mouse model of Alzheimer's disease[END_REF]. However, we shall see that the same is not true of models of molecular transport [START_REF] Goldman | Computational modeling of oxygen transport from complex capillary networks. Relation to the microcirculation physiome[END_REF][START_REF] Secomb | Green's function methods for analysis of oxygen delivery to tissue by microvascular networks[END_REF][START_REF] Fang | Oxygen advection and diffusion in a three-dimensional vascular anatomical network[END_REF][START_REF] Linninger | Cerebral microcirculation and oxygen tension in the human secondary cortex[END_REF][START_REF] Kojic | A composite smeared finite element for mass transport in capillary systems and biological tissue[END_REF][START_REF] Sweeney | Insights into cerebral haemodynamics and oxygenation utilising in vivo mural cell imaging and mathematical modelling[END_REF], which need to be improved. Most of them indeed rely on the well-mixed hypothesis and assume that radial gradients inside the blood vessels are negligible. The validity of this hypothesis requires detailed examination.

To illustrate this, let us first consider effective models of scalar transport in tubes, i.e. a larger class of problems that is central to many applications involving heat transfer, such as heat pipes in heat exchangers [START_REF] Faghri | Heat Pipe Science and Technology[END_REF]Reay, Kew & McGlen 2014) or solute transport, such as microchannels in labs-on-chip [START_REF] Bello | Use of Taylor-Aris dispersion for measurement of a solute diffusion coefficient in thin capillaries[END_REF][START_REF] Weigl | Microfluidic diffusion-based separation and detection[END_REF][START_REF] Beard | Taylor dispersion of a solute in a microfluidic channel[END_REF]. Since the influential works of [START_REF] Taylor | Dispersion of soluble matter in solvent flowing slowly through a tube[END_REF] and [START_REF] Aris | On the dispersion of a solute in a fluid flowing through a tube[END_REF], who focused on Poiseuille flow in impermeable straight cylinders, this problem has been extended to a large variety of configurations, including other velocity fields, e.g. [START_REF] Fan | Dispersion of Ostwald-de Waele fluid in laminar flow through a cylindrical tube[END_REF] and [START_REF] Gentile | The transport of nanoparticles in blood vessels: the effect of vessel permeability and blood rheology[END_REF], and more complex geometries, e.g. [START_REF] Brenner | Dispersion resulting from flow through spatially periodic porous media[END_REF] and [START_REF] Koch | Dispersion in fixed beds[END_REF]. Fundamental to scalar transport in a tube is the impact of radial gradients on the concentration field. These gradients, even if small in magnitude in the asymptotic regime, yield a huge increase in the longitudinal dispersion coefficient at sufficiently large Péclet numbers. For impermeable tubes, this is a consequence of the radial velocity gradient and the dispersion is often said to be shear induced or shear augmented. When the tube is permeable or semi-permeable, couplings between the scalar fields inside and outside the tube may yield additional gradients. In particular, sources or sinks in the domain surrounding the tube may generate strong radial concentration gradients within the tube. To study transport in the microcirculation, many authors [START_REF] Reneau | A digital simulation of transient oxygen transport in capillary-tissue systems (cerebral grey matter). Development of a numerical method for solution of transport equations describing coupled convection-diffusion systems[END_REF][START_REF] Levitt | Capillary-tissue exchange kinetics: an analysis of the Krogh cylinder model[END_REF][START_REF] Bate | Influence of diffusion on dispersion of indicators in blood flow[END_REF][START_REF] Lane | Indicator dilution measurement of mean transit time and flow in a straight tube[END_REF][START_REF] Hellums | The resistance to oxygen transport in the capillaries relative to that in the surrounding tissue[END_REF][START_REF] Tepper | Transient convective mass transfer in Krogh tissue cylinders[END_REF][START_REF] Lincoff | Characterisation of the unsteady transport of labelled species in permeable capillaries: role of convective dispersion[END_REF][START_REF] Baxter | Pharmacokinetic analysis of the perivascular distribution of bifunctional antibodies and haptens: comparison with experimental data[END_REF][START_REF] Fallon | Dispersion in core-annular flow with a solid annulus[END_REF][START_REF] Grinberg | Axial oxygen diffusion in the Krogh model[END_REF][START_REF] Secomb | Krogh-cylinder and infinite-domain models for washout of an inert diffusible solute from tissue[END_REF]) have considered a model axisymmetric system, called the Krogh cylinder, consisting of a single straight capillary surrounded by an annulus of tissue [START_REF] Krogh | The number and distribution of capillaries in muscles with calculations of the oxygen pressure head necessary for supplying the tissue[END_REF]. For this model system, results in [START_REF] Levitt | Capillary-tissue exchange kinetics: an analysis of the Krogh cylinder model[END_REF] and [START_REF] Lincoff | Characterisation of the unsteady transport of labelled species in permeable capillaries: role of convective dispersion[END_REF] suggest that Péclet numbers are small in the capillary bed and therefore that dispersive effects are negligible. However, the Krogh cylinder overly simplifies the geometry of the microvasculature so that extrapolation to realistic cases is difficult. For instance, cumulative effects may occur within networks, where dispersion and intravascular gradients might become important only for a large system. Further, since the Péclet number scales with both the characteristic length of the problem and the average velocity within each vessel, dispersive effects may become important in arterioles and venules [START_REF] Lincoff | Characterisation of the unsteady transport of labelled species in permeable capillaries: role of convective dispersion[END_REF]) that are larger in diameter and where blood flows faster.

We argue that recent developments in microscopy techniques are a key element in assessing the validity of the well-mixed hypothesis. For example, in vivo data in Sakadžić et al. (2014) show radial gradients of oxygen within arterioles, therefore experimentally demonstrating that the well-mixed hypothesis is not always valid. Other experiments also allow us to re-evaluate the importance of dispersive effects. Multiphoton microscopy [START_REF] Denk | Two-photon laser scanning fluorescence microscopy[END_REF][START_REF] Shih | Two-photon microscopy as a tool to study blood flow and neurovascular coupling in the rodent brain[END_REF], for example, can be used to determine the velocity within the brains of living rodents [START_REF] Santisakultarm | In vivo two-photon excited fluorescence microscopy reveals cardiac-and respiration-dependent pulsatile blood flow in cortical blood vessels in mice[END_REF][START_REF] Taylor | Microvascular basis for growth of small infarcts following occlusion of single penetrating arterioles in mouse cortex[END_REF] and infer the distribution of Péclet numbers in microvascular networks. Another approach consists in combining large post mortem anatomical network geometries [START_REF] Tsai | Correlations of neuronal and microvascular densities in murine cortex revealed by direct counting and colocalization of nuclei and vessels[END_REF]) with simulations of blood flow [START_REF] Schmid | Depth-dependent flow and pressure characteristics in cortical microvascular networks[END_REF][START_REF] Bracko | Neutrophil adhesion in brain capillaries contributes to cortical blood flow decreases and impaired memory function in a mouse model of Alzheimer's disease[END_REF] to calculate the velocity field throughout the network. Figure 2(b) shows the corresponding distribution of Péclet numbers for a molecule with high diffusivity (10 -9 m 2 s -1 ), which provides a lower bound for the Péclet numbers. It shows that most Péclet numbers are above 1 and suggests that, contrary to estimations made in [START_REF] Levitt | Capillary-tissue exchange kinetics: an analysis of the Krogh cylinder model[END_REF] and [START_REF] Lincoff | Characterisation of the unsteady transport of labelled species in permeable capillaries: role of convective dispersion[END_REF], Taylor's dispersion may be important in a significant proportion of the smallest capillaries.

In this work, we use a theoretical approach to ask whether the well-mixed hypothesis is adapted to describing molecular transport in the brain microcirculation. Our goal is to understand how radial concentration gradients affect intravascular transport, mass fluxes throughout the blood-brain barrier and molecular uptake within the tissue. Our strategy is to focus on situations of weak couplings between the vessels and tissue, so that the vascular network can be treated independently from the tissue ( § 2). We first use multiscale asymptotics to identify such situations. Then, we derive a generic form of the transport model for cross-section-averaged concentrations in tubes ( § 3.1) and networks ( § 3.2). We finally study the impact of dimensionless parameters and velocity profile on effective coefficients ( § 4.1) and compare our model with the well-mixed hypothesis at vessel ( § 4.2) and network scale ( § 4.3).

Problem set-up

Here, we set up the modelling bases of the problem of solute transport. Our goal is to derive a network model for complex microvascular architectures. To this end, we fi rst consider the case of a rigid, straight, long and narrow cylinder embedded in an infinite domain of tissue, as illustrated in fi gure 3. The solute can be metabolized within the tissue but behaves as a passive tracer in vessels, e.g. a nutrient. Following Pries et al. (1990) and [START_REF] Secomb | Blood flow in the microcirculation[END_REF], the blood is considered as a continuous fluid. The vessel walls are supposed to act as a membrane, semi-permeable to solute and impermeable to fl uid (Lorthois et al. 2014a;[START_REF] Kutuzov | Contributions of the glycocalyx, endothelium, and extravascular compartment to the blood-brain barrier[END_REF]. We also assume that the tissue surrounding the vessel is a continuous medium with uniform properties [START_REF] Levitt | Capillary-tissue exchange kinetics: an analysis of the Krogh cylinder model[END_REF] [START_REF] Nicholson | Diffusion and related transport mechanisms in brain tissue[END_REF][START_REF] Holter | Interstitial solute transport in 3D reconstructed neuropil occurs by diffusion rather than bulk flow[END_REF][START_REF] Kutuzov | Contributions of the glycocalyx, endothelium, and extravascular compartment to the blood-brain barrier[END_REF]. Since the aspect ratio (radius over length) of the vessel is small, the flow is treated as invariant along the vessel length and the velocity profi le is independent of the axial position along the vessel. Using these hypotheses, we fi rst write the mass balance equations within the vessel and tissue ( § § 2.1 and 2.2). Then we use multiscale asymptotics to describe the different transport regimes and identify situations of weak coupling where the vessel can be treated separately from the tissue ( § 2.3).

Microscopie transpon equations and boundary conditions

In order to derive the three-dimensional (3-D) local transport equations for the concentration within both vessels and tissue, we consider a binary solvent/solute mixture in each domain. The diffusion flux in vessels and tissue is determined using Fick's fi rst law. With these assumptions, the mass balance equation for solute within vessels can be written as Transport in the tissue is mainly diffusive so that we have

∂ t (ρ T Y T ) = ∇ • (ρ T D T ∇Y T ) -M, (2.2)
where ρ T is the local density of the tissue (kg m -3 ), Y T the local solute mass fraction, D T is the molecular diffusion coefficient (m 2 s -1 ) of the solute within the tissue and M the local metabolic reaction rate (kg m -3 s -1 ). Here, D V = D T and both are assumed to be uniform and constant. We consider that the solute is dilute and does not modify the densities. We can therefore write

C V = ρ V Y V , (2.3) C T = ρ T Y T , (2.4)
where ρ V and ρ T are uniform and constant and C V and C T are the local mass concentration (kg m -3 ). Hence, equations (2.1) and (2.2) now yield

∂ t C V = ∇ • (-UC V + D V ∇C V ), (2.5) ∂ t C T = ∇ • (D T ∇C T ) -M, (2.6)
which are the local transport equation for the concentration of solute within the vessels and tissue. The metabolic reaction rate M is modelled using a Michaelis-Menten kinetics, which is standard for oxygen [START_REF] Secomb | Green's function methods for analysis of oxygen delivery to tissue by microvascular networks[END_REF] or glucose [START_REF] Atkins | Physical Chemistry for the Life Sciences[END_REF]. Consequently, equation (2.6) becomes

∂ t C T = ∇ • (D T ∇C T ) -M T C T C 0 + C T , (2.7)
where M T is the maximal metabolic rate (kg m -3 s -1 ) and C 0 the concentration (kg m -3 ) at which the metabolic reaction term equals half its maximal value. Finally, the vessels and tissue are coupled together through the blood-brain barrier. We treat this barrier as a geometric manifold (infinitely thin layer) with transmission boundary conditions connecting the two domains. Assuming that there is no sorption of solute molecule onto the barrier, mass conservation reads

-n • (D V ∇C V ) = -n • (D T ∇C T ), (2.8)
where n is the unit normal vector pointing outward from the vessel. To represent the selective action of the blood-brain barrier, we use a membrane condition in the form

-n • (D V ∇C V ) = K m (C V -λC T ), (2.9)
with K m the membrane permeability (m s -1 ). Here, K m is primarily introduced as a parameter allowing us to model the ability of a given molecule to pass the barrier.

The limit K m → 0 corresponds to a molecule that cannot pass the barrier, resulting in a homogeneous Neumann boundary condition. On the other hand, K m → ∞ leads to C V = λC T , which corresponds to thermodynamic equilibrium and may be obtained through equality of chemical potentials with λ the partition coefficient [START_REF] Atkins | Physical Chemistry for the Life Sciences[END_REF]. This boundary condition, which is analogous to the thermal resistance in heat transfer, is widely used in microcirculation [START_REF] Lincoff | Characterisation of the unsteady transport of labelled species in permeable capillaries: role of convective dispersion[END_REF][START_REF] Fang | Oxygen advection and diffusion in a three-dimensional vascular anatomical network[END_REF][START_REF] Vikhansky | Taylor dispersion in finite-length capillaries[END_REF][START_REF] Sweeney | Insights into cerebral haemodynamics and oxygenation utilising in vivo mural cell imaging and mathematical modelling[END_REF]) and captures the selectivity of blood/tissue exchanges through the blood-brain barrier.

Non-dimensionalization of transport equations

The local transport problem, equations (2.5) and (2.7)-(2.9), is non-dimensionalized using the characteristic vessel length, L, and the associated diffusive time scale, L 2 /D V . The concentration field is non-dimensionalized using C inlet in the vessel and C 0 in the tissue. Here, C inlet corresponds to the concentration at the vessel inlet and C 0 to the concentration at which the metabolic reaction term equals half its maximal value. The dimensionless initial-boundary value problem now reads

∂ t C V = -Pe∇ • (U * C V ) + ∇ 2 C V ,
(2.10)

∂ t C T = η∇ 2 C T -ηDa T C T 1 + C T , (2.11) -n • ∇C V = -n • (ηγ ∇C T ),
(2.12) is the aspect ratio of the tube.

-n • ∇C V = Da m (C V -λγ C T ), (
2.3. Asymptotic transport regimes Due to the strong couplings between the vessel and tissue, it is not possible to derive a general 1-D effective transport equation for the cross-section-averaged concentration inside the vessel by simply averaging equation (2.10). As mentioned in the Introduction, one way to obtain such an effective model is to consider a subset of possible geometries, e.g. the Krogh cylinder. The main drawback of this approach is that it makes generalization to realistic microvascular architectures difficult. Another alternative is to identify regimes where couplings between the vessel and the tissue are relatively weak, so that the transport within each vessel can be treated as independent of the transport within the tissue. To this end, we use multiscale asymptotics [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Davit | Homogenization via formal multiscale asymptotics and volume averaging: how do the two techniques compare?[END_REF] to identify the specific scalings of the control parameters, expressed as powers of (e.g. Pe = O( i )), that yield weak couplings. In the following, we will further assume that, for small molecules like water, oxygen or nutrients, λ = O(1) and η = O(1). In addition, we suppose that the characteristic concentration in the tissue is small compared to the characteristic concentration inside the vessel, so that γ = O( ), as it is the case for oxygen [START_REF] Sweeney | Insights into cerebral haemodynamics and oxygenation utilising in vivo mural cell imaging and mathematical modelling[END_REF]. Doing so reduces the set of dimensionless numbers, for a given geometry, to Pe, Da T and Da m .

For the cylindrical geometry at hand (figure 3), the spatial differential operators have the following expressions:

∇• = e r ∂ r • +e z ∂ z •,
(2.16)

∇ 2 • = 1 r ∂ r (r∂ r •) + ∂ 2 z •, (2.17)
with e r and e z the unit vectors associated with the transverse and longitudinal directions. In such a configuration, there are primarily two scales of variations for the differential operators, the longitudinal scale z and the transverse scale r. In order to highlight those scales in (2.16) and (2.17), we rescale the radial coordinate and introduce the following change of variable:

ρ = r.
(2.18) Equations (2.16) and (2.17) then become

∇• = -1 e r ∂ ρ • +e z ∂ z •, (2.19) ∇ 2 • = -2 1 ρ ∂ ρ (ρ∂ ρ •) + ∂ 2 z • . (2.20)
For the rest of this section, we therefore define the following notations:

∇• = ∇ z • + -1 ∇ ρ • (2.21) and ∇ 2 • = ∇ 2 z • + -2 ∇ 2 ρ • . (2.22)
The main idea is then to separate the length scales of variation of the concentration fields in both vessel and tissue into different components and search for solutions in the form of series of powers of

C X = i i C X,i (r, z, t) X = V, T.
(2.23)

Injecting the above derivative expressions and the series form of the solution into the transport equations within vessel and tissues ((2.10) and (2.13)) yields

i i ∂ t C V,i = i i (-Pe(∇ z • (U * C V,i ) + -1 ∇ ρ • (U * C V,i )) + ∇ 2 z C V,i + -2 ∇ 2 ρ C V,i ), (2.24) i i ∂ t C T,i = i i     η∇ 2 z C T,i + -2 η∇ 2 ρ C T,i -ηDa T C T,i 1 + k k C T,k     .
(2.25) Doing the same for the boundary conditions at the vessel wall leads to

- i i n • (∇ z C V,i + -1 ∇ ρ C V,i ) = - i i n • ηγ (∇ z C T,i + -1 ∇ ρ C T,i ), (2.26) - i i n • (∇ z C V,i + -1 ∇ ρ C V,i ) = i i (Da m (C V,i -λγ C T,i )).
(2.27)

In (2.24), we have -1 ∇ ρ • (U * C V,i ) = 0 since there is only axial convection. Similarly, in (2.26) and (2.27) the terms n • ∇ z C V,i and n • ∇ z C T,i vanish as n (here e r ) is orthogonal to ∇ z • (here e z ∂ z •).

For each set of scalings of the dimensionless parameters (Pe, Da T , Da m ), we obtain a nested sequence of subproblems corresponding to the different powers of , which allows us to identify the leading-order concentration fields (C V,0 , C V,1 , C T,0 , C T,1 ). Effective transport equations for the cross-section-averaged concentration are then obtained by averaging the resulting equations in space via the following operator:

• = 1 π 2 •r dr dθ = 1 π •ρ dρ dθ (2.28) so that in general C V = i i C V,i . (2.29)
For concision, this cross-section-averaged concentration C V is simply referred to as the average concentration in the rest of this paper. Inspired by the work of [START_REF] Mei | Method of homogenization applied to dispersion in porous media[END_REF], [START_REF] Auriault | Taylor dispersion in porous media: analysis by multiple scale expansions[END_REF] and [START_REF] Allaire | Two-scale expansion with drift approach to the Taylor dispersion for reactive transport through porous media[END_REF], we illustrate how this approach can be used to obtain the well-known regime of Taylor's dispersion. We show that, in this case, the boundary conditions ((2.12) and (2.13)) degenerate into Neuman boundary conditions, consistent with the impermeable walls of Taylor's case ( § 2.3.1). Then, we go on to generalize this study to a large range of parameter scalings and to identify the regimes in phase space that correspond to weak vessel-tissue couplings ( § 2.3.2). For all these scalings, we further show that the membrane boundary condition (2.13) degenerates into a Robin-type boundary condition at the vessel walls.

Recovering Taylor's dispersion

Taylor's dispersion is obtained when molecules primarily remain inside the vessels and the Péclet number is large enough. The corresponding set of scalings reads

Pe = O( -1 ), Da m = O( 2 ), Da T = O( -3 ).
(2.30a-c)

Because the Péclet number is large, it is necessary to decompose the time derivative within the vessel using two time scales to capture both the contributions of longitudinal convection and longitudinal diffusion to the effective transport. This reads

∂ t • = ∂ t * • +Pe∂ τ •, (2.31)
where t * corresponds to the time variation of longitudinal diffusion and τ corresponds to the time variation of longitudinal convection.

In addition, Da T = O( -3 ) implies that the tissue Damköhler is very large. In such a regime, we expect the concentration in the tissue to be small so that the reaction rate can be approximated by a first-order kinetics. Note that taking into account the nonlinearity in the reaction rate does not change the results and only leads to tedious additional calculations that are not necessary for understanding the procedure. For this set of parameters, in the limit → 0, equations (2.24) and (2.25) lead to

C T,0 = 0 O( -3 ),
(2.32) (2.37) and equation (2.27)

∇ 2 ρ C V,0 = 0 η∇ 2 ρ C T,0 -η( 3 Da T )C T,1 = 0 O( -2 ), (2.33) ( Pe)∂ τ C V,0 = -( Pe)∇ z • (U * C V,0 ) + ∇ 2 ρ C V,1 η∇ 2 ρ C T,1 -η( 3 Da T )C T,2 = 0 O( -1 ), (2.34) ∂ t * C V,0 + ( Pe)∂ τ C V,1 = -( Pe)∇ z • (U * C V,1 ) + ∇ 2 z C V,0 + ∇ 2 ρ C V,2 ∂ t C T,0 = η∇ 2 z C T,0 + η∇ 2 ρ C T,2 -η( 3 Da T )C T,3 O(1). (2.35) Equation (2.26) becomes -n • ∇ ρ C V,0 = 0, (2.36) -n • ∇ ρ C V,i = -η( -1 γ )n • ∇ ρ C T,i-1 , i > 0,
-n • ∇ ρ C V,i = 0, i 2, (2.38) -n • ∇ ρ C V,3 = ( -2 Da m )C V,0 , i = 3.
(2.39)

-n • ∇ ρ C V,i = ( -2 Da m )(C V,i-3 -λ( -1 γ )C T,i-4 ), i > 3.
(2.40)

By mathematical induction, we obtain that C T,i = 0 ∀i ∈ N inside the tissue, so that the problem is fully controlled by intravascular transport. Equation (2.38) further shows that, at leading order, there is no diffusive flux at the vessel/tissue interface so that the membrane boundary condition degenerates into a homogeneous Neumann condition. Therefore, the model is equivalent to an impermeable wall, which obviously removes all vessel/tissue couplings and corresponds to Taylor's dispersion regime. The effective transport equation is obtained by averaging equation (2.35) using the operator defined in (2.28). This leads to (2.41) which can be further simplified by noticing that ∇ 2 ρ C V,2 = 0 using the divergence theorem and the boundary condition at the vessel wall (2.38). To obtain an expression for U * C V,1 in (2.41), we consider equation (2.33) and use the divergence theorem to obtain C V,0 = C V,0 (z). Therefore, equation (2.34) can be written as

∂ t * C V,0 + Pe∂ τ C V,1 = -Pe∇ z • U * C V,1 + ∇ 2 z C V,0 + ∇ 2 ρ C V,2 ,
( Pe)∂ τ C V,0 = -( Pe)U * (ρ) • ∇ z C V,0 + ∇ 2 ρ C V,1 . (2.42)
Averaging equation (2.34) leads to

( Pe)∂ τ C V,0 = -( Pe) U * • ∇ z C V,0 . (2.43)
We have already shown that

C V,0 = C V,0 (z), which implies that C V,0 = C V,0 .
Consequently, we can subtract (2.43) from (2.42) to obtain

( Pe) Ũ(ρ) • ∇ z C V,0 = ∇ 2 ρ C V,1 , (2.44)
where Ũ = U * -U * is the local perturbation of the velocity field. We then have a solution in the form

C V,1 = ( Pe)b(ρ) • ∇ z C V,0 , (2.45)
where b is the closure variable which solves

Ũ = ∇ 2 ρ b, (2.46) n • ∇ ρ b = 0, (2.47)
Modelling solute transport in the brain microcirculation with uniqueness obtained via the average constraint b = 0.

(2.48)

In this regime, the definition of the average concentration implies that

C V = C V,0 , C V,1 = 0, (2.49)
so that injecting (2.45) into (2.41) leads to

∂ t * C V = (1 -2 Pe 2 U * • b )∇ 2 z C V . (2.50) Now, substituting ∂ t * C V using (2.31) yields ∂ t C V = Pe∂ τ C V + (1 -2 Pe 2 U * • b )∇ 2 z C V .
(2.51)

Finally, using (2.43) in the above expression leads to

∂ t C V = -Pe U * • ∇ z C V + (1 -2 Pe 2 U * • b )∇ 2 z C V , (2.52)
which is the effective transport equation, where (1 -2 Pe 2 U * • b ) represents an effective diffusion coefficient for Taylor's regime with a general expression of the velocity field U * . In the remainder of this work, we will call this the generalized Taylor's dispersion model. The closure variable b captures the radial gradients of concentration that stem from gradients in the velocity profile. The expression -2 Pe 2 U *

• b in (2.52) characterizes hydrodynamical dispersion effects. As expected, the dispersion coefficient scales as the square of the Péclet number, potentially leading to differences of several orders of magnitude in the spreading of average concentration when compared to the well-mixed hypothesis (C V ≈ C V ). The latter is equivalent to C V,1 = 0, yielding b = 0, and U *

• b = 0 in (2.52), and therefore cancelling all dispersive effects.

Phase diagram: identification of weak couplings

The approach illustrated for the generalized Taylor's dispersion in the above section has been repeated for a large range of scalings of Pe, Da m , Da T . Although the developments are tedious, they are relatively systematic and well described in the literature (e.g. see [START_REF] Koch | Dispersion in fixed beds[END_REF][START_REF] Mei | Method of homogenization applied to dispersion in porous media[END_REF][START_REF] Auriault | Taylor dispersion in porous media: analysis by multiple scale expansions[END_REF][START_REF] Auriault | Transport in porous media: upscaling by multiscale asymptotic expansions[END_REF]. Therefore, we only present here a summary of the results obtained, in the form of phase diagrams of transport regimes. The key point is that, for some specific scalings, the membrane boundary condition can be simplified into a classic Robin boundary condition, which we will show corresponds to weak vessel/tissue couplings. These scalings are represented by the dotted regions in the phase diagram of transport regimes, displayed in the parameter plane (Pe, Da m ) for Da T = O( -3 ) in figure 4(a) and in the parameter plane (Da T , Da m ) for Pe = O(1) in figure 4(b). They include a variety of cases that will be described thereafter, such as impermeable cases or situations where the concentration within the tissue is very small because the reaction in the tissue is almost instantaneous. Among these regions of weak couplings, some scalings lead to equations that can be homogenized using multiscale asymptotics, and the others to equations that cannot be homogenized. An example of such a region is the case where Pe

= O( -2 ), Da m = O( 2 ) and Da T = O( -3 ). For these ´4 ´-5 ´3 ´2 ´1 1 Pe ´-1 ´-2 ´-3 ´-4 ´-5 ´4 ´3 ´2 ´1 1 ´-1 ´-2 ´-3 ´-4 ´-5 ´-4 ´-3 ´-2 ´-1 Da m Da T 1 ´1 ´2 ´3 ´4
´-5

´-4

´-3

´-2

´ scalings, multiscale asymptotic analysis shows that both C V,0 and C V,1 are spatially homogeneous. The associated transport equation for the average concentration yields

∂ t C V,0 = -Pe∇ z • ( U * C V,2 ),
(2.53)

where U * C V,2 cannot be made explicit without writing all the remaining transport equations corresponding to the higher-order concentration fields (C V,1 , C V,2 ...), making it impossible to homogenize, as discussed in [START_REF] Auriault | Taylor dispersion in porous media: analysis by multiple scale expansions[END_REF]. Below, we thus focus on the homogenizable regions. These have been highlighted in green in figure 4. Each roman numeral corresponds to a different class of transport equation involving various effective coefficients, such as the effective diffusion coefficient introduced in (2.52). The expressions for these effective coefficients are given in appendix A so that, in the following, we only focus on the structure of transport equations.

For small values of the membrane Damköhler number (Da m = O( i ), i 2), the membrane boundary condition degenerates towards a Neumann condition at leading orders ( j 2)

n • ∇ ρ C V,j = 0, j i, (2.54) n • ∇ ρ C V,j = ( -i Da m )C V,j-(i+1) , j = i + 1. (2.55) n • ∇ ρ C V,j = ( -i Da m )(C V,j-(i+1) -λ( -1 γ )C T,j-(i+2) ), j > i + 1. (2.56)
Hence, in these regimes, intravascular transport is fully decoupled from the tissue and solely depends on the Péclet number. For a low Péclet number (Pe = O( i ), i 1), transport is purely driven by molecular diffusion (domain I in figure 4a) and convection plays a negligible role. In this regime the effective transport equation is Domain I :

∂ t C V = ∇ 2 z C V .
(2.57)

When the Péclet number increases (Pe = O(1)), convection and molecular diffusion play similar roles (domain III in figure 4a,b) so that the effective transport equation reads Domain III :

∂ t C V = -Pe∇ z C V + ∇ 2 z C V .
(2.58)

For even larger values of the Péclet number (Pe = O( -1 )), shear-enhanced diffusion dominates (generalized Taylor's regime, domain V in figure 4a) and the influence of radial gradients of concentration cannot be neglected (see § 2.3.1). The effective transport equation is

Domain V : ∂ t C V = -Pe∇ z C V + D eff (Pe)∇ 2 z C V , (2.59)
where the effective diffusion coefficient D eff (Pe) is strictly larger than one. For Pe = O( i ), i < -1, it is no longer possible to formally derive an effective transport equation for the average concentration, as discussed before for the case Pe = O( -2 ).

For larger membrane Damköhler numbers (Da m = O( i ), i < 2), the vessel walls are no longer impermeable. Starting with Da m = O( ), the membrane boundary condition degenerates towards Neumann and Robin conditions at leading orders ( j 2)

n • ∇ ρ C V,j = 0, j 1, (2.60) n • ∇ ρ C V,2 = ( -1 Da m )C V,0 , j = 2. (2.61) n • ∇ ρ C V,j = ( -1 Da m )(C V,j-(i+1) -λ( -1 γ )C T,j-(i+2) ), j > 2,
(2.62) so that the transport in the vessel remains independent of the transport inside the tissue. Then, for Da m = O( i ), i < 1, the transport in the vessel becomes strongly coupled to the tissue and, in general, the problem cannot be restricted to intravascular transport. The only exception is when the Damköhler number in the tissue is large (Da T = O( i ), i -3). In such regimes, all molecules passing through the vessel walls are instantaneously metabolized, which always results in C T being negligible (similar to § 2.3.1). In this situation, the membrane condition reads

n • ∇ ρ C V,j = 0, j i, (2.63) n • ∇ ρ C V,j = ( -i Da m )C V,j-(i+1) , j > i, (2.64) in the range Da m = O( i ), -1 i < 1 and C V,j = 0, j -(i + 2), (2.65) n • ∇ ρ C V,j+i+1 = ( -i Da m )C V,j , j > -(i + 2), (2.66) in the range Da m = O( i ), i < -1.
= O( i ) -1 < i < 2, domain II in figure 4a) yields Domain II : ∂ t C V = ∇ 2 z C V -K eff (Da m ) C V , (2.67)
where K eff (Da m ) > 0 is the effective reaction rate representing the exchanges with the tissues. To obtain (2.67), for the specific scaling Da m = O(1), it is necessary to introduce an additional time scale so that 4a,b) means that we have

∂ t • = ∂ t * • + -1 Da m ∂ ν
= O( i ) -1 < i < 2, domain IV in figure
Domain IV : ∂ t C V = -Pe∇ z C V + ∇ 2 z C V -K eff (Da m ) C V .
(2.68)

Again, to obtain (2.68), for the specific scaling Da m = O(1), it is necessary to use a time scale decomposition analogous to the one introduced for domain II.

For larger Péclet numbers (Pe = O( -1 ), domain V), increasing the membrane Damköhler number (Da m = O( ), domain VI in figure 4a) also results in the appearance of an effective reaction rate

Domain VI : ∂ t C V = -Pe∇ z C V + D eff (Pe)∇ 2 z C V -K eff (Da m ) C V , (2.69)
and increasing further the membrane Damköhler number (Da m = O(1), domain VII in figure 4a) introduces an additional effective convection term,

Domain VII : ∂ t C V = -U eff (Pe, Da m )∇ z C V + D eff (Pe)∇ 2 z C V -K eff (Da m ) C V ,
(2.70) where U eff (Pe, Da m ) > 0 is the effective velocity. In order to obtain (2.70), it is necessary to introduce two additional time scales so that

∂ t • = ∂ t * • +Pe∂ τ • + -1 Da m ∂ ν • to
account for both longitudinal convection and radial diffusion through the blood-brain barrier.

For very large Péclet numbers (Pe = O( i ), i < -1), it is again not possible to derive an effective transport equation for the average concentration. Similarly, it is not possible to derive such an equation for a very large membrane Damköhler number (Da m = O( i ), i -1).

All the above results have been obtained by assuming that the reference concentration in the tissue was small compared to the reference concentration in the vessel so that γ = O( ). For the case where the reference concentrations are comparable, i.e. when γ = O(1), we find similar results except that, in the (Da T , Da m ) plane, domain IV is only valid for Da m = O( ) when strong reaction takes place in the tissue, i.e. Da T = O( i ), i -3.

In summary, in all cases of weak couplings (figure 4, dotted areas), the membrane boundary condition can be simplified into a Robin condition that captures the Neumann case for Da m = O( i ), i 2 ((2.54) and (2.55)), the Robin case for Da m = O( i ), -1 i < 2 ((2.63) and (2.64)) and the Dirichlet case for Da m = O( i ), i < -1 ((2.65) and (2.66)). The corresponding local transport equations are therefore equivalent to

∂ t C V = -Pe∇ • (UC V ) + ∇ 2 C V ,
(2.71)

-n • ∇C V = Da m C V .
(2.72)

These equations can only be homogenized by multiscale asymptotics for the specific scalings highlighted in green in figure 4, and correspond to different classes of effective transport equations (Domains I-VII).

Altogether, equations (2.71) and (2.72) model situations for which the transport problem inside the vasculature is independent of the concentration field in the tissue and the concentration field in the tissue is controlled by the concentration field within the vessel. In such situations, the concentration in the tissue can be reconstructed a posteriori, for instance using Green's functions approaches (see [START_REF] Hsu | A Green's function method for analysis of oxygen delivery to tissue by microvascular networks[END_REF], based on the knowledge of the flux at the vessel wall ((2.12) and (2.72)). In that sense, vessels and tissue are weakly coupled.

Effective transport model for weak couplings

While the existence of various classes of transport equations is not an issue for effective models at the scale of a single vessel, it becomes problematic when studying transport in microvascular networks. Indeed, a wide range of parameter scalings co-exist in networks, as illustrated in figure 2(b) for Péclet numbers, and using the equations above would imply switching between the different classes of equations for each vessel.

To overcome this difficulty, the main idea of this section is that the regimes of validity of the homogenizable Robin boundary condition having been identified in figure 4 (green regions), we now can use (2.71) and (2.72) as a starting point to derive a unified effective transport equation encompassing all the above classes (Domains I-VII). Getting inspiration from the more general problem of solute transport in porous media with surface reaction, e.g. [START_REF] Shapiro | Taylor dispersion of chemically reactive species: irreversible first-order reactions in bulk and on boundaries[END_REF], [START_REF] Golfier | Heat and mass transfer in tubes: an analysis using the method of volume averaging[END_REF] and [START_REF] Van Duijn | Chapter 1. Effective dispersion equations for reactive flows with dominant Péclet and Damkohler numbers[END_REF], we adopt the method of volume averaging with closure [START_REF] Whitaker | The Method of Volume Averaging[END_REF][START_REF] Davit | Homogenization via formal multiscale asymptotics and volume averaging: how do the two techniques compare?[END_REF]) for a single vessel ( § 3.1) and generalize it to microvascular networks ( § 3.2).

Unified transport equation for single vessels

Spatial averaging

From now on, we will use the radial coordinate r instead of ρ with r = ρ to describe the transport in the transverse direction. Therefore, now only represents the dimensionless radius of the vessel, so that r = corresponds to the vessel wall and C V ( , z, t) is the concentration at the wall. The differential operators ∇• and ∇ 2

• are applied following the definitions given in (2.16) and (2.17). With these definitions, we can proceed to the first step of volume averaging, which consists in averaging equation (2.71) using the operator defined in (2.28). This leads to

∂ t C V = -Pe∇ • U * C V + ∇ 2 C V + 1 π (n • ∇C V ) dθ , (3.1)
where (1/π ) (n • ∇C V ) dθ is the cross-section-averaged mass flux across the bloodbrain barrier stemming from averaging the radial component of the Laplacian operator.

The above equation constitutes the basis for the unified effective transport equation with only U * C V and (1/π ) (n • ∇C V ) dθ left to be made explicit as functions of C V . To do so, we use a perturbation decomposition for both the concentration and velocity fields,

C V = C V (z) + CV (r, z), (3.2) U * = U * + Ũ(r), (3.3)
where the velocity field perturbation Ũ is the same as the velocity field introduced in (2.44). Applying the average operator on (3.2) and (3.3) yields

CV = 0, (3.4) Ũ = 0, (3.5)
since U * and C V do not depend on radial position, and averaging them does not change their values. Introducing the perturbation decomposition for the concentration field into (2.71) and (2.72) yields

∂ t C V + ∂ t CV = -Pe∇ • (U * C V ) -Pe∇ • (U * CV ) + ∇ 2 C V + ∇ 2 CV , (3.6) -n • ∇ C V -n • ∇ CV = Da m C V + Da m CV . (3.7)
Subtracting (3.1) from (3.6) leads to

∂ t CV = -Pe∇ • ( Ũ C V + U * CV -Ũ CV ) + ∇ 2 CV - 1 π n • (∇ C V + ∇ CV ) dθ , (3.8) -n • ∇ C V -n • ∇ CV = Da m C V + Da m CV . (3.9)
Assuming that the perturbation field relaxes much faster than the average fields (quasistationarity, see discussions in [START_REF] Davit | Homogenization via formal multiscale asymptotics and volume averaging: how do the two techniques compare?[END_REF]), we can write

-Pe∇ • ( Ũ C V + U * CV -Ũ CV ) + ∇ 2 CV + 1 π Da m ( C V + CV ) dθ = 0. (3.10)
At this stage, to solve the solute transport problem within the vessel, one would have to solve both the average transport problem (3.1) and the perturbation problem formed by (3.9) and (3.10). To uncouple these equations, a closure must be introduced.

The closure problem

The following order-one closure is used:

CV = α(r) C V (z) + β(r)∂ z C V (z), (3.11)
with α and β two functions of the radial position, r. Since variables depending on r and z are now separated, it is possible to write the perturbation problem in cylindrical coordinates. Substituting the closure in (3.9) and (3.10) and further assuming that second-order derivatives can be neglected leads to

(-Pe Ũ -PeU * α + Pe Ũα + β + 2 -1 Da m β( ))∂ z C V + (α + 2 -1 Da m (1 + α( ))) C V = 0 (3.12) and (α ( ) + Da m (1 + α( ))) C V + (β ( ) + Da m β( ))∂ z C V = 0, ( 3
∂ t C V = -U eff ∂ z C V + D eff ∂ 2 z C V -K eff C V , (3.22)
with the effective velocity

U eff = Pe + Pe αU * + 2 -1 Da m β( ), (3.23)
the effective diffusion coefficient

D eff = 1 -PeU * β (3.24)
and the effective reaction rate

K eff = 2 -1 Da m (1 + α( )). (3.25)
We note that D eff has an expression similar to the effective diffusion coefficient introduced in (2.52) for the general Taylor's dispersion regime, except that here PeU * β describes the effect of radial gradient of concentration on axial dispersion induced by both the velocity profile and the molecular exchanges. Consequently, equations (3.22)-(3.25) can be seen as the unified version of the effective transport equations for the average concentration. In the following, these equations will be referred to as the weak coupling averaged (WCA) model. The effective coefficients in this model depend on the Péclet and the membrane Damköhler numbers. They also depend on the shape of the velocity profile, which controls the interactions between the velocity and concentration gradients at microscopic scale.

For comparison, we are now deriving the effective transport equation using the well-mixed hypothesis, for which we neglect radial gradients of concentration, so that C V ≈ C V . Similarly to the volume averaging procedure, we start by averaging equation (2.71) which leads to

∂ t C V = -Pe∇ • U * C V + ∇ 2 C V + 1 π (n • ∇C V ) dθ , (3.26)
where n = e r is the unit vector associated with radial direction. Now, using the decomposition of the concentration field introduced in (3.2) for the last term of the above equation leads to

∂ t C V = -Pe∇ • U * C V + ∇ 2 C V + 1 π (n • ∇ C V + n • ∇ CV ) dθ . (3.27)
By definition of the averaging operator, we know that

C V = C V (z). Therefore we have n ⊥ ∇ C V , leading to n • ∇ C V = 0.
In addition, since we assume that C V ≈ C V , we can deduce that U * C V ≈ U * C V , turning the above equation into

∂ t C V = -Pe∇ • ( U * C V ) + ∇ 2 C V + 1 π (n • ∇ CV ) dθ . (3.28)
The same analysis can be done starting from the Robin condition at the vessel wall (3.7), which leads to

-n • ∇ CV = Da m C V . (3.29)
Finally, equations (3.28) and (3.29) can be combined to yield

∂ t C V = -Pe∂ z C V + ∂ 2 z C V -2 -1 Da m C V . (3.30)
This is the well-mixed version of the effective transport equation for the average concentration, for which we have In these limits, we have α = 0 and β = 0 within the vessel, α = 0 and β = 0 at the vessel wall, and therefore α = β = 0 with the average condition. This leads to

U WM = Pe, (3.31) D WM = 1, (3.32) K WM = 2 -1 Da m . ( 3
U eff -U WM = Pe αU * + 2 -1 Da m β( ), (3.34) D eff -D WM = -PeU * β , (3.35) K eff -K WM = 2 -1 Da m α( ), (3.36)
which means that the well-mixed hypothesis not only leads to misestimating of the dispersion coefficient, but also the other effective parameters.

Extension to microvascular networks

Here, we extend the WCA and WM models from a single tube to a network. In § 3.2.1, we start by quickly presenting the blood flow model that allows us to predict the velocity in each vessel of a microvascular network. Then, in § 3.2.2, we use these velocity fields to derive the constitutive equations of molecular transport for the particular case of a bifurcation and show that they can be used for larger network geometries.

Blood flow

The vasculature is treated as a network of interconnected tubes, where the stationary flow rate distribution is computed using a nonlinear network approach described in [START_REF] Pries | Blood flow in microvascular networks. Experiments and simulation[END_REF], [START_REF] Lorthois | Simulation study of brain blood flow regulation by intra-cortical arterioles in an anatomically accurate large human vascular network. Part I. Methodology and baseline flow[END_REF] and [START_REF] Bracko | Neutrophil adhesion in brain capillaries contributes to cortical blood flow decreases and impaired memory function in a mouse model of Alzheimer's disease[END_REF]. This approach, where blood is considered as a homogeneous fluid and red blood cells are treated as a volume fraction (haematocrit), accounts for the complex rheological properties of blood flow in microcirculation through two in vivo empirical laws. The first one, which accounts for the Fåhraeus-Lindquist effect, describes the average dissipation at vessel scale through an apparent viscosity which depends on the tube diameter and haematocrit, so that a linear relationship between the flow rate and the pressure drop can be written in each tube [START_REF] Pries | Microvascular blood viscosity in vivo and the endothelial surface layer[END_REF]. The distribution of haematocrit in the network and phase-separation effects are captured by the second empirical law [START_REF] Pries | Red cell distribution at microvascular bifurcations[END_REF]) that links haematocrit and flow rate ratios at bifurcations. This problem is nonlinear and is solved iteratively (see e.g. [START_REF] Lorthois | Simulation study of brain blood flow regulation by intra-cortical arterioles in an anatomically accurate large human vascular network. Part I. Methodology and baseline flow[END_REF][START_REF] Fry | Estimation of blood flow rates in large microvascular networks[END_REF][START_REF] Bracko | Neutrophil adhesion in brain capillaries contributes to cortical blood flow decreases and impaired memory function in a mouse model of Alzheimer's disease[END_REF] for given boundary conditions (e.g. imposed pressures at network inlets and outlets, imposed haematocrit at network inlets), yielding the pressure, the flow rate and the haematocrit within the network.

Molecular transport

To extend the effective transport problems derived in § 3.1.3 to microvascular networks, the relationships between the average concentration at the ends of connected tubes must be modelled. In the brain, these networks are mainly composed of bifurcations connecting three vessels. Thus, we only detail here the specific case of a single bifurcation with three vessels, schematized in figure 2(a). We start by writing the effective transport equation for the average concentration within each vessel of the bifurcation

(∂ t C V + U eff ∂ z C V -D eff ∂ 2 z C V + K eff C V ) i = 0, (3.37)
where i ∈ {0, 1, 2} denotes the vessel index for the three vessels connected to a central bifurcation vertex, B, which we consider to be a material point with no volume, consistent with the 1-D form of the WCA effective equation. At the bifurcation vertex, the average concentration is assumed to be continuous. Therefore, we introduce C B the average concentration at the bifurcation (vertex B) and C V i (z B ) the concentration in vessel i in the limit z → z B , the position of B. We have

C V i (z B ) = C B i ∈ {0, 1, 2}. (3.38)
The total molecular flux conservation is also imposed as

i n i • (U app C V e z -D eff ∂ z C V e z ) i (z B ) dS i = 0, i ∈ {0, 1, 2}. (3.39)
Here n i and dS i are the normal vector and elementary surface associated with the cross-section of vessel i, and

U app = U eff -2 -1 Da m β( ) (3.40)
is the apparent velocity, which is different from the effective velocity. The term 2 -1 Da m β( ) in the effective velocity (3.23) indeed directly comes from the exchange term in (3.1), hence it models molecules exiting the vessel through the blood-brain barrier. However, when considering the intravascular flux, only the flow velocity corresponding to convection must be used. This can be proven by considering the physical expression of the flow in terms of C V , not C V , and using the closure relationships. Finally, to close the system formed by (3.37)-(3.39), we impose standard boundary conditions, such as a Dirichlet or Neumann conditions at the inlets and outlets. This approach, illustrated for a simple bifurcation, can be extended to an entire network in a straightforward manner.

Results

Our goal is to apply the above developments to complex networks. For that purpose, we first need to evaluate the range of dimensionless parameters encountered in the brain microcirculation. Values of Pe can be estimated from figure 2(b), where the radial Péclet number ( Pe) is typically between 10 -3 and 10 2 for small, highly diffusible vital molecules such as oxygen. These numbers increase by an order of magnitude for a solute with a tenfold lower diffusivity, as is the case for gadolinium [START_REF] Benson | Threedimensional concentration field measurements in a mixing layer using magnetic resonance imaging[END_REF][START_REF] Wieseotte | An estimate of Gd-DOTA diffusivity in blood by direct NMR diffusion measurement of its hydrodynamic analogue Ga-DOTA[END_REF], an intravascular contrast agent frequently used in clinical imaging. The radial Damköhler number ( Da m ) typically varies between 0 for such an intravascular tracer to an upper bound value corresponding to small, vital molecules diffusing easily through the blood-brain barrier. For such molecules, an upper bound for the membrane permeability can be roughly estimated as K m = D V /e, where e is the thickness of the endothelial layer ( 10 -6 m, [START_REF] Kutuzov | Contributions of the glycocalyx, endothelium, and extravascular compartment to the blood-brain barrier[END_REF]. This yields Da m = R/e, i.e. larger than 1 in all vessels. For this range of parameters, we study the behaviour of U eff , D eff and K eff as a function of Pe and Da m ( § 4.1) for each vessel separately. Together, these parameters control the average concentration in the WCA model. To further assess their impact, we compare this model with the well-mixed model at vessel scale ( § 4.2) and network scale ( § 4.3).

Impact of the dimensionless parameters and velocity profile on effective

coefficients Here, we derive explicit expressions of U eff , D eff and K eff for a generic polynomial expression of the velocity profile ( § 4.1) and study in detail the case of a twoparameter velocity profile including both bluntness and slip velocity ( § 4.1.2), which are characteristic of microvascular blood flows (see e.g. [START_REF] Damiano | Estimation of viscosity profiles using velocimetry data from parallel flows of linearly viscous fluids: application to microvascular haemodynamics[END_REF][START_REF] Roman | Velocimetry of red blood cells in microvessels by the dual-slit method: effect of velocity gradients[END_REF][START_REF] Lei | Blood flow in small tubes: quantifying the transition to the non-continuum regime[END_REF][START_REF] Roman | Going beyond 20 µm-sized channels for studying red blood cell phase separation in microfluidic bifurcations[END_REF].

General polynomial shape

We consider the following polynomial approximation for the velocity profile:

U * = i=N i=0 w i r i , (4.1)
with w i chosen so that U * = 1. This expression allows us to solve the closure problems (3.16)-(3.21) and obtain analytical expressions of α and β (see appendix B). Using (3.23)-(3.25), we then obtain analytical expressions of U eff , D eff and K eff as functions of the aspect ratio and the Péclet and membrane Damköhler numbers.

The effective velocity varies linearly with the Péclet number and can be written as

U eff = Pe(1 + U + ), (4.2)
where U + is an apparent overspeed defined by

U + = - 4 Da m Da m + 4 i w i i (i + 2)(i + 4) + 2( Da m ) 2 ( Da m + 4) 2 i w i i 2 + 2i + 8 (i + 2)(i + 4)(i + 6)
.

(4.3) We will show later that U + is positive for all velocity fields considered in § 4.1.2, a posteriori justifying its denomination. This overspeed only depends on the membrane Damköhler number and the aspect ratio (4.3). Recalling that Da m , the radial membrane Damköhler number, represents the magnitude of vessel/tissue exchanges compared to diffusion at microscopic scale suggests that the apparent overspeed is directly linked to radial gradients of concentration resulting from the vessel/tissue exchanges. Consistent with this idea, we have Da m = 0 for an impermeable tube yielding U + = 0 for any velocity profile. As a result, the effective velocity is equal to the Péclet number, as expected from the phase diagram in figure 4. Moreover, in the diffusion-limited regime, when Da m → ∞, U + always remains bounded, which is critical for consistency at large membrane Damköhler numbers.

The effective diffusion coefficient can be written as

D eff = 1 + ( Pe) 2 Pe 2 c , (4.4)
where Pe is the radial Péclet number and Pe c is a critical Péclet number, strictly positive and bounded for all velocity fields considered later, defined by

Pe -2 c = -2 i j w i w j (i + 2)( j + 2) j (i + 4)(i + j + 4) - j 4( j + 4) - 2 Da m Da m + 4 i j w i w j (i + 2)( j + 2)( j + 4) 2( j + 2) (i + j + 6) - j + 4 (i + j + 4) - 2j (i + 4) + ( j + 2) 2 2( j + 6) + 2( Da m ) 2 ( Da m + 4) 2 i j w i w j ( j 2 + 2j + 8)i 4(i + 2)(i + 4)( j + 2)( j + 4)( j + 6) . (4.5)
For a given velocity profile, the critical Péclet number only depends on the radial membrane Damköhler number. As expected, in the case of Taylor's dispersion ( Da m = 0, w 0 = 2, w 1 = 0, w 2 = -2, N = 2), the critical Péclet number is equal to √ 48.

The critical Péclet number can be interpreted in two different ways. First, it represents a regime threshold in terms of radial Péclet number. When Pe Pe c , the contribution of radial gradients of concentration to the axial effective diffusion becomes larger than the contribution of molecular diffusion. Alternatively, it can be viewed as the sum of the contributions of different effects. Equation (4.4) can indeed be rewritten as

D eff = 1 + ( Pe) 2 Λ impermeable + ( Pe) 2 Λ couplings ( Da m ), (4.6)
where Λ impermeable is independent of the membrane Damköhler number and represents the contribution of the radial gradients in the case of an impermeable membrane, while Λ couplings captures the coupling between the velocity profile and exchanges at the vessel walls.

Finally, the effective reaction rate is

K eff = 8 -1 Da m Da m + 4 , (4.7)
which does not depend on the velocity profile and monotonically increases with the membrane Damköhler number. Once again, this effective reaction rate is bounded when Da m → ∞, which is critical for consistency in the diffusion-limited regime.

Two-parameter velocity profile

We now consider the following two-parameter velocity profile

U * r = (n + 2) k(n + 2) -2 k - r n , (4.8)
with k 1, n 2. The parameter k controls the slip velocity and the parameter n controls the bluntness of the velocity profile. This form of velocity profile allows us to capture both the slip (red blood cells can roll along the walls, see [START_REF] Roman | Velocimetry of red blood cells in microvessels by the dual-slit method: effect of velocity gradients[END_REF][START_REF] Sherwood | Spatial distributions of red blood cells significantly alter local haemodynamics[END_REF]) and the perturbation of the velocity profile induced by the presence of the red blood cells (flatter profiles may result from the migration of red blood cells towards the centre of the vessel, see e.g. [START_REF] Pries | Biophysical aspects of blood flow in the microvasculature[END_REF][START_REF] Santisakultarm | In vivo two-photon excited fluorescence microscopy reveals cardiac-and respiration-dependent pulsatile blood flow in cortical blood vessels in mice[END_REF][START_REF] Roman | Going beyond 20 µm-sized channels for studying red blood cell phase separation in microfluidic bifurcations[END_REF].

Figure 5 shows the velocity as a function of the radial position for different values of (k, n). The particular case k = 1, n = 2 corresponds to Poiseuille flow. Increasing either the slip velocity or the profile bluntness results in flattening the velocity profile near the vessel centreline, making it more uniform. Using such a parametrization allows us to investigate in more detail the physical behaviour of the apparent overspeed and of the critical Péclet number.

Figure 6 displays U + as a function of the radial membrane Damköhler number for different values of (k, n). Increasing the membrane Damköhler number for a given (k, n) monotonically increases the apparent overspeed, from zero (impermeable walls) to an upper-bound value which depends on both k and n (diffusion-limited regime). Slower molecules close to the vessel walls indeed have a higher chance of leaving the bloodstream and entering the surrounding tissue. On the contrary, for a given Damköhler number, increasing the slip velocity or the velocity profile bluntness results in the velocity being more uniform across the vessel cross-section, which reduces the relative importance of low speed areas close to the walls and decreases the apparent overspeed. Thus, the largest apparent overspeed is obtained for Poiseuille flow in the limit of large Da m . In this case, U + → 1, meaning that the contribution (k,n). Because the critical Péclet number captures an interplay between two different contributions (see (4.6)), its variations are non-monotonie with a minimum that depends on the shape of the velocity profile and two asymptotes in the limits EDa 111 --* 0 and EDa 111 --* oo. In the Iimit EDa 111 --* 0 (impermeable tube), Pe c is constant and the effective diffusion is driven by Taylor's dispersion with Pe c = .J4s for Poiseuille flow. In the diffusion-limited regime corresponding to EDa 111 --* oo, the concentration vanishes at the vesse] wall (see § 4.2.1), so that the low velocities close to the vesse] walls do not contribute as much to shear-induced dispersion. In other words, the molecules that remain in the tube effectively explore a Jess heterogeneous velocity fi eld and the critical Péclet number is increased. For a given Damkohler number, increasing the slip velocity or the velocity profi le bluntness increases the critical Péclet number.

1 5 ••••••••••• -••••• •• • •• .. ..

5�-----------�

10-3 10-2 10-1 Hl" 10 1 102 103

EDa ,. Figures 6 and7 show that the effective velocity and effective diffusion coefficient are smaller for (k > 1, n > 2) than for a Poiseuille profi le. However, equations ( 4.2) and (4.4) imply that they are larger than the WM ones for any velocity profile. In the same way, equation (4.7) implies that the effective reaction rate is always smaller than the WM rate. Therefore, as soon as the radial membrane Damkohler number reaches 1, radial gradients play a signifi cant role on transport at the scale of a single vesse!, even when the overspeed and effective axial diffusion are relatively small due to the shape of the velocity profile. The impact on both transport and blood/tissue exchanges at the scale of a single vesse! will be studied in the next section.

Comparison of the WCA model with well-mixed models at vessel scale

Here, we investigate how these differences in the effective coefficients impact both transport and blood/tissue exchanges at the scale of a single vesse!. To this end, we compare the WM and WCA models for the case of Poiseuille flow in both steady and transient regimes. Corresponding solutions for the two models are derived analytically in appendix C. To further validate the WCA model, 2-D axisymmetric numerical solutions of the local transport problem ((2.71) and (2.72)) are obtained via finite volume (see appendix D for details).

Stationary regime

In this section, the inlet concentration is always set to one and the outlet effective diffusive flux set to zero. Figure 8 displays the average concentration along the vesse! axis for different radial Péclet and membrane Damkohler numbers, with insets showing the associated radial concentration profile. The fi rst important result is that the WCA model (green continuous line) is in excellent agreement with the reference numerical solutions (blue dash-dotted line) for a wide range of Péclet and membrane Damkohler numbers. We see that the WM model is also in excellent agreement with this reference solution when the radial membrane Damkohler number is small (fi gure 8a,c). For Pe = 0 (fi gure 8a), this is because both shear-induced and exchange effects are negligible. For high Péclet numbers (fi gure 8c), this is because the average concentration along the vesse! axis is alrnost constant as a result of our choice of inlet/outlet boundary conditions and the low exchanges with the tissues. In this specific case, the radial gradients of concentration are not significant The situation is fundamentally different for large Damkohler numbers. In this case, the well-mixed model underestimates the average concentration and overestimates the decrease along the vessel (figure 8b,d). Exchanges with the tissue are strong and radial concentration gradients are important (see figure 8b,d insets). These gradients are not accounted for in the well-mixed model, so that the flux from the blood to the tissue is overestimated. Increasing the Péclet number further amplifies this issue because of additional shear-induced gradients (figure 8d).

We now study how these differences impact blood/tissue exchanges. To quantify these, we calculate the integral flux through the vessel wall normalized by the inlet mass flux cross the blood-brain barrier (E = 1) is obtained for rather small values of the radial membrane Damkohler number (�10-1 ) (fi gure 9a). This is because, despite the low membrane permeability, all molecules have enough time to cross the blood-brain barrier before being convected out of the vesse!. Conversely, for large Péclet numbers, the WM mode! signifi cantly overestimates the extraction coefficient for moderate to large Damkohler numbers (fi gure 9b). In these regimes, convection is so strong that a signifi cant fraction of molecules reaches the vesse! outlet and exchanges become limited by radial diffusion. By contrast, the WM mode! assumes that molecular diffusion in the cross-section is instantaneous, which overestimates the effective reaction rate. The latter is even unbounded for large membrane Damkohler numbers, as shown by (3.33). As a consequence, for an arbitrarily large value of the inlet mass flux, the WM mode! always predicts a Damkohler number for which E = 1, which is physically inconsistent.

E = J J (Uapp(C) -DeffV (C)) • ndS;n1e1 -J J (Uapp(C) -D eff V (C)) • ndSoiùtet !! , (4.9) (Uapp(C) -D eff V (C)) • ndSinlet

Transient transpon of an initial square function

We now focus on an idealized transient regime where the tracer initially fills up 20 % of the vesse! length (black dotted line in fi gure 10) to represent solute injection. The temporal evolution of this initial square function is also displayed in figure 10 for two values of the Damkohler number. Consistent with the above stationary results, we find that the WCA mode! is always in excellent agreement with the reference numerical solution. For an impermeable tube (figure 10a, EDa m = 0), the differences between the WM and WCA predictions are only due to Ta ylor's dispersion. This explains the smaller spreading of the concentration pulse for the well-mixed case, which only considers molecular diffusion. For a permeable tube (fi gure 10b, EDa m = 0.5), the WM mode! predicts a concentration maximum upstream the reference solution. This is a direct consequence of the apparent overspeed, which is not captured by the WM model.

4.3.

Comparison of the WCA and well-mixed models in a large microvascular network In this section, our goal is to investigate the cumulative effects of radial concentration gradients at the scale of a large anatornical rnicrovascular network previously obtained from a mouse brain by [START_REF] Tsai | Correlations of neuronal and microvascular densities in murine cortex revealed by direct counting and colocalization of nuclei and vessels[END_REF] and [START_REF] Blinder | The cortical angiome: an interconnected vascular network with noncolumnar patterns of blood flow[END_REF]. To this end, we solve equations (3.37)-(3.39) in the network shown in figure 2(a), with effective coefficients corresponding either to the WCA mode] (with k = 1, n = 2) or the WM mode], in transient and stationary regimes ( § § 4.3.1 and 4.3.2, respectively). The flow rate distribution in this network is computed using the approach described in § 3.2.1, enabling us to deterrnine the velocity distribution (fi gure 2b) needed to solve the transport problem. This transport problem is further closed using the following boundary conditions: transient or stationary Dirichlet conditions at all arteriolar inlets; natural Neumann conditions at all venular outlets; pseudo-periodic boundary conditions on the lateral faces (see details in Cruz-Hernandez et al. ( 2019)); natural Neumann conditions on capillary vessels crossing the bottom face. Such conditions ensure that, for impermeable walls, ail molecules entering through the arteriolar inlets exit through the venular outlets. The whole transport problem is solved analytically (see appendix C) for steady regimes and numerically using fi nite volume (see appendix D) for transient regimes. ln large networks, such solutions are too complex to be interpreted by focusing on the behaviour of the average concentration along the vesse] axes, as illustrated by fi gure 11 and supplementary movie 1 (found at https://doi.org/10.1017/jfm.2019.866) in a transient regime. To overcome this issue, we thus only consider integral quantities at the scale of the network. Results are presented below for two classes of molecules: (i) small molecules with high diffusivity (Dv = 10-9 m 2 s-1 ), such as water, free ox ygen or glucose; and (ii) larger molecules with lower diffusivity (Dv = 10-w m 2 s-1 ), such as gadolinium, which is frequently used in clinical applications.

Transient regimes

We first study the transient injection of a tracer that does not cross the blood-brain barrier (intravascular tracer, Da m = 0) by considering square input functions of variable durations. We quantify the model response by calculating the transient integral output flux

I(t) = all outlets (U app C -D eff ∇ C ) • n dS. (4.10)
Figure 12 displays this flux as a function of time for the two classes of molecules and for different input durations, represented by bold lines on the top x-axis (0.01 s in figure 12a,b and 0.1 in figure 12c,d).

Regardless of the effective transport model, the molecular diffusion coefficient or the duration of the injection, we find that the integral output flux increases sharply after a delay corresponding to the smallest transit time between inlets and outlets throughout the network. After this delay, which is always longer for the well-mixed case, the behaviour of the output flux is strongly influenced by the effective transport model, except for long injection durations and the largest diffusivity (figure 12d). The WM model exhibits multiple peaks with fast temporal dynamics (red dotted lines), while the WCA model exhibits a much smoother behaviour (plain green lines) with a single local maximum. When Taylor's dispersion is neglected (WM model, red dotted lines), the multiple peaks are a signature of the multiple convective pathways connecting each arteriolar inlet to several venular outlets [START_REF] Lorthois | Simulation study of brain blood flow regulation by intra-cortical arterioles in an anatomically accurate large human vascular network. Part I. Methodology and baseline flow[END_REF][START_REF] Guibert | Coupling and robustness of intracortical vascular territories[END_REF]. Sorne of these pathways have short transit times (of order 0.1 s), either because they are short and/or are associated with high velocities (preferential pathways), while others have much longer transit times (of order 1 s). When the duration of the injection is smaller than the transit time difference between two pathways, these appear as two separate peaks, unless molecular diffusion is sufficiently strong to merge the two peaks. Thus, either increasing the injection duration (figure 12c) or increasing the diffusivity (fi gure 12b) results in decreasing the number of peaks. In the same way, accounting for Taylor's dispersion increases the effective diffusion coefficient and always significantly smooths out the signal (red dotted versus plain green lines in all panels). However, when Taylor's dispersion is accounted for, increasing the diffusion coefficient results in steeper variations, regardless of the injection duration. This counter-intuitive result is due to Taylor's dispersion depending on the square of the Péclet number. By increasing the diffusion coefficient, the Péclet number is decreased and so is the effective diffusion coefficient, therefore causing steeper variations in the integral output fl ux. Finally, we observe that the long-term relaxation of the output flux is not affected by Taylor's dispersion and that the tails of both the WM and WCA model are almost identical for times greater than approximately 0.5 s. This is because the tail corresponds to pathways associated with smaller Péclet numbers, where Taylor's dispersion becomes negligible. In these pathways, the concentration signal bas more time to spread due to molecular diffusion, which explains the gentle tailing behaviour. Similarly, for very long injection durations (5 s), the two models yield almost identical results, as shown in figure 13. To summarize, figure 12 shows that, for a purely intravascular tracer (Da m = 0), the influence of Taylor's dispersion is the strongest when molecular diffusion is small and injection duration is short.

Stationary regimes

We now focus on blood/tissue exchanges in a stationary regime with permeable vessel walls, similar to the regime studied for single vessels in § 4.2.1. To reduce the number of parameters, we fix the membrane permeability so that the radial Damkohler number is the same for all vessels, for any value of the aspect ratio. This might slightly affect the details of the results, as aspects ratio are known to depend on vessel type (e.g. venules slightly larger than arterioles, see Lorthois, Lauwers & Cassot 2014b). However, we have checked that this does not change the general conclusions presented below. The extraction coefficient is defined for arterioles, venules and capillaries as This allows us to distinguish the respective contribution of arterioles, venules and capillary vessels to blood/tissue exchange at network scale. The total network extraction is obtained by summing up these separate contributions.

Figure 14 shows that the total extraction coefficient increases with the radial membrane Damköhler number (black lines). In the asymptotic regime of large Damköhler numbers, the total extraction coefficient reaches a plateau slightly smaller than one for low diffusivity molecules and equal to one for highly diffusive molecules. In the latter case, all molecules cross the blood-brain barrier before reaching any network outlet. For both classes of molecules, the stationary extraction coefficient associated with the arterioles (red lines) monotonically increases with the Damköhler number. Arterioles are indeed the most upstream vessels of the network and are perfused with a constant input of solute concentration. Thus, increasing the membrane Damköhler number always increases the arterioles capacity to exchange molecules with the surrounding tissue. By contrast, capillary extraction behaves non-monotonically with the Damköhler number. For small Damköhler numbers, the capillary domain is the main contributor to the total extraction coefficient. In this regime, vessel walls are nearly impermeable and only a few molecules leave the vessels to enter the tissue (E tot 1). The whole network is perfused with a concentration that is almost constant and the total extraction coefficient is mainly controlled by the largest exchange surface, i.e. the capillary bed. For large Damköhler numbers, however, upstream vessels (arterioles) exchange more with the tissue and only transmit a fraction of the input mass to downstream vessels (capillaries and venules). As a consequence, the capillary bed is not perfused as well as in the small Damköhler number regime and the capillary extraction coefficient is smaller. This transition between good and partial perfusion explains the existence of a maximum for E cap in all cases. A similar analysis holds for the venules, the most downstream vessels, so that their associated extraction coefficient increases only for small Damköhler numbers where the network is well perfused. For large Damköhler numbers, the venule extraction coefficient either reaches a small value plateau because venules are perfused with a residual concentration, as shown by figure 14(a), top left, or goes to zero because arterioles and capillaries have already exchanged all molecules, as shown in figure 14(b).

Consistent with results presented in § 4.2.1, this global blood/tissue exchange behaviour is well captured by the WM model for small radial Damköhler numbers, regardless of the vessel territory or the molecular diffusion coefficient (figure 14c,d). In this regime, the exchanges are weak and the concentration is globally homogeneous throughout the network. For large Damköhler numbers, however, the WM model overestimates the intensity of exchanges (with E tot = 1 regardless of the molecular diffusion coefficient). As a result, the contribution of upstream vessels is strongly overestimated and arterioles are the sole contributor to the exchanges at the largest Damköhler numbers. In this model, as previously discussed, radial diffusion is instantaneous and the exchanges are never limited by molecular diffusion, ultimately causing all molecules to cross the blood-brain barrier before even exiting the arteriolar domain.

Discussion

Since the influential papers on solute dispersion in impermeable tubes by [START_REF] Taylor | Dispersion of soluble matter in solvent flowing slowly through a tube[END_REF] and [START_REF] Aris | On the dispersion of a solute in a fluid flowing through a tube[END_REF], the huge impact of vanishingly small radial concentration gradients generated by non-uniform velocity fields, and the associated physics, has been well understood. In these geometries, the resulting shear-induced dispersion usually scales with the square of the Péclet number, so that the effective diffusion term in the corresponding 1-D averaged equation is, in principle, neither negligible at small Péclet numbers (where molecular diffusion dominates) nor at large Péclet numbers (where dispersion becomes large compared to advection). However, due to the mathematical and computational difficulties in accounting for all transport phenomena at the scale of large microvascular networks, many of the recent papers modelling oxygen supply to the brain assume purely axial convection within the vessels [START_REF] Safaeian | A computational model of oxygen transport in the cerebrocapillary levels for normal and pathologic brain function[END_REF][START_REF] Gagnon | Quantifying the microvascular origin of BOLD-fMRI from first principles with two-photon microscopy and an oxygen-sensitive nanoprobe[END_REF][START_REF] Gould | The capillary bed offers the largest hemodynamic resistance to the cortical blood supply[END_REF][START_REF] Sweeney | Insights into cerebral haemodynamics and oxygenation utilising in vivo mural cell imaging and mathematical modelling[END_REF]) and use a well-mixed hypothesis for radial diffusion. This hypothesis neglects radial gradients of concentration altogether, whether they originate from heterogeneities in the velocity fields or other mechanisms. This has broad implications, as it eliminates the impact of gradients on all the effective parameters such as dispersion or vessel/tissue exchange coefficients. This idea of perfect mixing has been recently contradicted by experimental measurements that show significant radial oxygen gradients in the penetrating arterioles of living anesthetized mice (Sakadžić et al. 2014). This suggests the existence of couplings between the vessels and tissue, which are not captured in Taylor's problem due to the impermeable walls, but are likely to have a significant impact on solute transport within the brain vasculature. Such couplings are difficult to study in the general case of complex interconnected microvascular networks. In such networks, it is indeed challenging to define a geometry representative of the tissue in which it would be possible to derive effective equations describing transport in the tissue.

Here, our strategy was to take a step back in modelling this complex problem and to focus on asymptotic regimes where vessels and tissue are only weakly coupled. Using multiscale asymptotics, we demonstrated that these regimes of weak couplings, while involving various effective transport equations, are all associated with situations where the solute concentration is negligible at all times and at all locations within the tissue. This occurs either because the consumption rate is sufficiently strong and/or the exchange rate is sufficiently small. In these cases, we demonstrated that the boundary condition at the vessel walls degenerates into a classic Robin boundary condition. Of course, this situation is encountered only for a very limited subset of solutes with direct interest either in physiology or in medical imaging, including purely intravascular tracers (tracers that do not cross the blood-brain barrier). This restricts the direct applicability of the present model to specific problems encountered in real life. However, the mechanisms at play, in particular radial gradients generated by the velocity field and by outward fluxes, are also present when stronger couplings are at play. Thus, this model approach allows us to answer the fundamental question of whether radial gradients can be generally ignored or should be treated carefully in transport models of the brain microcirculation. To address this question and identify the corresponding regimes, we compared the results of a new effective model, obtained using volume averaging, and the model based on the well-mixed hypothesis. To further assess the accuracy of our model, we also systematically compared these results against fully resolved numerical solutions of the microscopic transport problem. The main results obtained by this strategy can be summarized and put in context as follows.

For transient regimes, we primarily focused on the case where the blood-brain barrier is impermeable to solutes and tracers are purely intravascular. We demonstrated that Taylor's dispersion has a strong impact on tracer dynamics at network scale for injections durations typically smaller than the longer transit time of blood from arterioles to venules, approximately 1 s for the network considered in this study. This may be important to improve the interpretation of perfusion imaging techniques that aim at quantifying blood flow in various organs, including the brain. The general idea of such techniques is to transiently perfuse a tracer, whether a radio-labelled inert gas, radio-labelled water or magnetized water, through the vasculature. In arterial-spin labelling magnetic resonance angiography (ASL-MRA), for example, the labelling duration used to magnetize water protons is typically approximately 1.8 s [START_REF] Alsop | Recommended implementation of arterial spin-labeled perfusion MRI for clinical applications: a consensus of the ISMRM perfusion study group and the European consortium for ASL in dementia[END_REF]. The dynamics of the transport is subsequently measured. An image of the distribution of blood flow is then reconstructed by solving an inverse problem based on a simplified direct model of transport within the brain. When the characteristic times are approximately 1 s or smaller, Taylor's dispersion will significantly affect the direct model, the resolution of the inverse problem, and therefore the image of blood flow.

Of course, this is just an order of magnitude as our model is limited and other important effects will need to be evaluated in the future. For instance, the blood-brain barrier is not impermeable to water protons. Therefore, could the apparent overspeed affect the distribution of transit times throughout the network? Also, we considered volumes of approximately 1 cubic millimetre, which is small in comparison with the total volume of the cortex (∼1000 cubic millimetres in rodents, Kovačević et al. 2005). Then, does the size of the volume modify the dynamics sufficiently so that dispersive effects needs to be accounted for in ASL-MRA? Finally, how do strong vessel/tissue couplings, for which Robin boundary conditions cannot be used, affect the characteristic times, especially for the longer injection durations used in other perfusion imaging techniques, where vessel/tissue equilibrium is expected (Lorthois et al. 2014a)?

For steady regimes, we considered a case with blood/tissue exchanges. This situation is important to address fundamental questions in physiology such as oxygen delivery and its interactions with brain function or with mechanisms of vascular development in health and disease, e.g. the hypoxia-driven development of vessels in tumours [START_REF] Leonard | The analysis of convection and diffusion in capillary beds[END_REF]. At network scale, we showed that radial concentration gradients strongly impact the magnitude of the exchanges when the radial Damköhler number is above one. This regime is the only one where significant concentration gradients are expected to build up within the vessels. Therefore, the fact that radial oxygen gradients measured experimentally are very strong suggests that the radial Damköhler number for oxygen must be above one, consistent with the rough estimates presented in § 4. This likely to be also the case for other small vital molecules such as water and glucose.

For both transient and steady regimes, we found that at the scale of single vessels, the well-mixed hypothesis should not be used when both the Péclet and the Damköhler numbers are above one. In cerebral microvascular networks, we already know that radial Péclet numbers are above one in many vessels for a wide range of solutes in physiological conditions. Therefore, the well-mixed hypothesis is not valid at network scale as soon as the radial Damköhler number is larger than one in a significant number of vessels.

Based on these results, we believe that it is important to explore in more detail the role of spatio-temporal heterogeneities on transport through vascular systems. Indeed, many other types of heterogeneities have been neglected. For instance, the above results have been obtained for a homogeneous fluid but blood is not homogeneous. It is a dense suspension of red blood cells with a complex microvascular flow structuration (e.g. single line red blood cell flow in capillary vessels, cell free layer). One of the most widely characterized consequence of such a structuration is a deviation from the Poiseuille flow profile in small micro-vessels with diameters below ∼100 µm. In this paper, we described such effects using a time-averaged velocity profile that includes bluntness and slip [START_REF] Damiano | Estimation of viscosity profiles using velocimetry data from parallel flows of linearly viscous fluids: application to microvascular haemodynamics[END_REF][START_REF] Roman | Velocimetry of red blood cells in microvessels by the dual-slit method: effect of velocity gradients[END_REF][START_REF] Lei | Blood flow in small tubes: quantifying the transition to the non-continuum regime[END_REF][START_REF] Roman | Going beyond 20 µm-sized channels for studying red blood cell phase separation in microfluidic bifurcations[END_REF]. We demonstrated that both phenomena contribute to decreasing the impact of radial concentration gradients on microvascular solute transport. However, we know very little about the large-scale implications of temporal heterogeneities and of considering a time-average velocity profile. Further, the influence of other gradient-generating mechanisms should also be considered. For example, restricted diffusion of specific tracers through red blood cell membranes could decrease the effective radial diffusion coefficient in blood and therefore increase the gradient magnitude. On the other hand, shear-induced radial dispersion of red blood cells [START_REF] Bishop | Effect of aggregation and shear rate on the dispersion of red blood cells flowing in venules[END_REF][START_REF] Kabacaoglu | Quantification of mixing in vesicle suspensions using numerical simulations in two dimensions[END_REF][START_REF] Tang | Shear-induced diffusion of red blood cells measured with dynamic light scattering-optical coherence tomography[END_REF] may have the opposite effect. How, then, does that affect transport at the scale of the brain?

Conclusion

In this work, we have presented a new effective model of solute transport in large microvascular networks. The novelty of this model is that it does not neglect radial gradients of concentration, but rather is valid for all regimes where the boundary condition at vessel walls can be assimilated to a Robin condition (regimes of weak vessel/tissue couplings). Our simulations show that models based upon the well-mixed hypothesis fail to accurately capture the spreading of solute concentration in transient problems for processes faster than approximately 1 Hertz. In stationary regimes, these models also wrongly evaluate the amount and distribution of solute mass exchanged between vessels and tissue for molecules that easily cross the blood-brain barrier. This demonstrates that radial gradients are a fundamental component of transport in such systems and opens the way towards a more global characterization of the impact of spatio-temporal heterogeneities on transport in the microcirculation. This may ultimately lead the way towards a better understanding of transport in physiological and pathological conditions, or even developing more accurate models for the interpretation of perfusion imaging techniques. C.1. Single vessel Altogether, equations (3.22) and (3.30) are simple 1-D time-dependent equations, which can be solved analytically in a single vessel, provided that the effective coefficients are known and that the vessel inlet/outlet boundary conditions are simple enough. Starting with the stationary regime and assuming that D eff > 0 and K eff > 0 , the general solution of (3.22) can be written as

r i+2 -1 - 1 2(i + 2) r 2 -1 - i 4(i + 2)(i +
C V (z) = Ae ω 1 z + Be ω 2 z , ( C 1)
with ω 1 and ω 2 defined as follows:

ω 1 = U eff -U 2 eff + 4D eff K eff 2D eff ω 2 = U eff + U 2 eff + 4D eff K eff 2D eff , (C 2a,b)
and A and B coefficients to be determined from inlet and outlet boundary conditions. For example, for an inlet concentration equal to one (Dirichlet condition) and an effective diffusive flux at the outlet equal to zero (homogeneous Neumann condition)

C V stat (z) = ω 1 ω 2 e ω 1 +ω 2 ω 2 e ω 2 -ω 1 e ω 1
e ω 1 (z-1) ω 1 - where w 1 and u>i are defined in (C 2) and the coefficients A; and B; are to be determined using the boundary conditions at the inlet and outlet of each vessel. To do so, the above expression is substituted into (3.38) and (3.39), which gives (Ae "' • ZB + Be wzZB ); = (C) 8 , i E {O, 1, 2}, for the diffusive flux have been chosen. To validate this numerical framework, the stationary analytical solution introduced in appendix A (C 6)-(C 8) is obtained in the network presented in figure 2(a), for different molecule diffusivities D V and membrane permeabilities K m . As an example, figure 15 displays the errors between the finite volume resolution and this analytical solution as a function of the spatial discretization (degrees of freedom), for D V = 1e -9 m 2 s -1 and K m determined in each vessel so that Da m = 0.1. Among the three error estimators the L 2 error is the slowest to converge with roughly a slope of -1, since the advection transport scheme (QUICK), initially of order 3, degenerates towards a classic order-one upwind scheme for degrees of freedom neighbouring bifurcation vertices. Such a result validates the numerical framework described above.

e ω 2 (z-1) ω 2 . ( C

FIGURE 1 .

 1 FIGURE 1. Schematics of molecular transport processes at different scales in the brain with (a) scale of the whole brain; (b) scale of the cortex; (c) scale of a microvascular network; (d) scale of a single capillary vesse!. This fi gure is a modified version of fi gure 1 in Lorthois et al. (2014a).

FIGURE 3 .

 3 FIGURE 3. Schematic of solute transport within a single capillary. Here Cv is the local concentration of solute within the vessel and Cr the local concentration within the tissue. U(r) represents the velocity profile; Dv the diffusion coefficient of the solute within the vessel; K m is the membrane permeability; Dr denotes the diffusion coefficient within the tissue; and M (Cr) the consumption rate of solute within the tissue.

  (2.1) where Pv is the local mass density in the vessel (kg m-3 ), Yv the local solute mass fraction, U the local velocity (m s-1 ) and Dv the molecular diffusion coefficient (m 2 s-1 ) of the solute in vessels.

  2.13) where Pe = U L/D V is the longitudinal Péclet number based on the cross-sectionaveraged velocity U = (1/πR 2 ) U • n dS and U * = U/ U is the normalized velocity; Da T = M T L 2 /C 0 D T is the longitudinal Damköhler number within the tissue; η = D T /D V is the diffusivity ratio; Da m = K m L/D V is the longitudinal membrane Damköhler number; and γ = C 0 /C inlet is the ratio of reference concentrations. The corresponding radial Péclet and Damköhler numbers, based on the cylinder radius instead of its length, are Pe radial = Pe, Da T,radial = 2 Da T , Da m,radial = Da m ,

  FIGURE 4. Phase diagrams of transport regimes in (a) the (Pe, Da m ) plane for Da T = O( -3 ) and in (b) the (Da T , Da m ) plane for Pe = O(1). Dotted areas correspond to regimes of weak couplings where the membrane condition becomes a classic Robin condition. Green areas, which highlight the main result of this figure, correspond to parameter regimes for which it is possible to derive an effective transport equation for the cross-section-averaged concentration by multi-scale asymptotics. The associated transport regimes are labelled by roman numerals with: I, molecular diffusion only; II, molecular diffusion with effective reaction; III, advection and molecular diffusion; IV, advection, molecular diffusion with effective reaction; V, generalized Taylor's dispersion; VI, advection, effective diffusion and effective reaction; VII, effective advection, effective diffusion and effective reaction.

  12) and (3.13) must remain valid for any value of C V or ∂ z C V . The closure variables must thus solve the following set of ordinary differential equations:α + 2 -1 Da m (1 + α( )) = 0, (3.16) -Pe Ũ -PeU * α + Pe Ũα + β + 2 -1 Da m β( )formed by (3.16)-(3.21) defines the closure problem, which can be solved analytically for specific velocity profiles (e.g. a generic polynomial profile see § 4.1 and appendix B). The closure variables are then injected into equation (3.1) in order to make explicit U * C V and (1/π ) (n • ∇C V ) dθ . 3.1.3. The effective transport equations (weak coupling averaged and well-mixed models) Injecting the solution of the closure problem into the averaged transport equation (3.1) leads to

  30)-(3.33) will be referred to as the well-mixed (WM) model. These expressions are consistent with previous developments as they can be directly obtained by solving the closure problem defined by (3.16)-(3.21) when Pe → 0 and Da m → 0.

FIGURE 5 .

 5 FIGURE 5. Blood velocity profile as a function of (a) the parameter n that controls the profile bluntness (for no slip velocity, k = 1), and (b) the parameter k that controls the slip velocity (for a parabolic profile, n = 2).

FIGURE 6 .

 6 FIGURE 6. Overspeed as a function of the radial membrane Damkohler number for different (a) velocity profi le bluntness and (b) slip velocities.

  . -. .. . -' ...... . .. . .. . .. . -• .. . ... .. .... -.. ....

  .................... .. 10 ... .. .. . ..... .... .. , ... .. ,--• . .. . ..

FIGURE 7 .

 7 FIGURE 7. Critical Péclet number as a function of radial membrane Damkohler number for different (a) velocity profile bluntness and (b) slip velocities.

FIGURE 8 .

 8 FIGURE 8. Average stationary axial concentration for different Péclet and membrane Damkohler numbers for E=Ü.05. (a) EPe = O and EDa m = O.l. (b) EPe = O and EDa m = lO. (c) EPe = 100 and EDa m = 0.1. (d) EPe = 100 and EDa m = 10. Green lines are analytical solutions of the WCA model. Red dashed lines are analytical solutions of the WM model. Blue dash-dotted lines are full 2-D finite volume solutions of the local transport problem lnsets show radial concentration profiles.

FIGURE 9 .

 9 FIGURE 9. Extraction coefficient as a function of the radial membrane Damkohler number for E = 0.05. (a) EPe = 1. (b) EPe = 100. Green lines are analytical solutions of the WCA model. Red dashed lines are analytical solutions of the WM model. Blue dash-dotted lines are full 2-D finite volume solutions of the local transport problem

FIGURE 10 .

 10 FIGURE 10. Transport of an initial square function at different times for E = 0.1. (a) EPe = 5 and EDa m = 0. (b) EPe= 5 and EDa m = 0.5. Green lines are analytical solutions of the WCA model. Red dashed lines are analytical solutions of the WM model. Blue dash-dotted line s are full 2-D fi nite volume solutions of the local transport problem .

FIGURE 11 .

 11 FIGURE 11. Numerical solutions of average concentration of an intravascular tracer (Da m = 0) with a high diffusion coefficient (D V = 10 -9 m 2 s -1 ) resulting from a square input of 0.1 s duration at different times: (a) T = 0.1 s; (b) T = 0.2 s; (c) T = 0.3 s; (d) T = 0.5 s; (e) T = 0.6 s; ( f ) T = 0.7 s. This solution has been obtained using the WCA model on the network shown in figure 2(a) and using the velocity field computed in Cruz-Hernández et al. (2019). The case presented here corresponds to the plain green curves presented in figure 12(d). For full transient dynamics (T = 0 s to T = 1 s), see supplementary movie 1.

FIGURE 12 .

 12 FIGURE 12. Integral output mass flux as a function of time: (a) 0.01 second input Dv = 10-10 m 2 s-1 ; (b) 0.01 second input Dv = 10-9 m 2 s-1 ; (c) 0.1 second input Dv = 10-10 m 2 S-1 ; {a) 0.1 second input Dv = 10-9 m 2 S-1 • Green lines are analytical solutions of the WCA model. Red dashed lines are analytical solutions of the WM model. Thick black lines represent the input duration.

FIGURE 13 .

 13 FIGURE 13. Integral output mass fl ux for a 5 s input square function Dv = 10-9 m 2 s-1 . Sirnilar results have been obtained for Dv = 10-10 m 2 s-1 • The green line corresponds to the WCA model, the red dashed line to the WM model. The thick black line represents the input duration.

�FIGURE 14 .

 14 FIGURE 14. Stationary network extraction coefficient as a function of radial membrane Damköhler number. (a) WCA model and D V = 10 -10 m 2 s -1 . (b) WCA model and D V = 10 -9 m 2 s -1 . (c) WM model and D V = 10 -10 m 2 s -1 . (d) WM model and D V = 10 -9 m 2 s -1 . Tot., total extraction coefficient; Art., arteriolar extraction coefficient; Ven., venular extraction coefficient; Cap., capillary extraction coefficient.

  Both stationary and transient analytical solutions of the effective transport equation in a single vessel (3.22) are first given for classical inlet and outlet boundary conditions. Similarly, stationary solutions of the associated system of equations at the scale of the network (3.37)-(3.39) are then derived.

FIGURE 15 .

 15 FIGURE 15. Error estimators for the comparison between reference analytical and numerical solutions of (3.37}-(3.39) in the large anatomical anatomie network displayed in figure 2(a), as a function of the number of degrees of freedom (DoF) for Dv = 10-9 m 2 s-1 and EDam = 0.1. Error estimators: emean = (1/NooF) I: ICvF -c,efl, emax = maxooF ICvF -C,etl, eL2 = IICvF -C,etll-

(

  new closed linear problem for A; and B;, provided that boundary conditions at the inlets and outlet of the bifurcation are known. Once A; and B; are determined for each vesse], the concentration throughout the bifurcation can be deduced using expression (C 6). Appendix D. Numerical methods for solving the effective transport equation ln transient regimes, equations (3.37} -(3.39) are solved numerically using fülite volume. Briefly, given the simple expression of the effective transport equation and the range of dimensionless parameter involved in microcirculation, a simple classic QUICK (order-3 upwind) scheme for the convective flux and a order-2 centred scheme

et aL 2017; Sweeney et al. 2018

; Lincoff et al. 1983; Hellums et al. 1995; Goldman & Popel 1999; Secomb et al. 2004; Fallon & Anuj 2005; Fang et al. 2008; Linninger et al. 2013; Kojic ) where transport is mainly diffusive

  •, where ∂ ν represents the time scale associated with radial diffusion through the blood-brain barrier. For a moderate Péclet number (Pe = O(1), domain III), increasing the membrane Damköhler number (Da m

  Appendix C. Analytical solutions of the effective transport equations

	+ 2 Pe	Da m Da m + 4		i	w i		i (i + 2)(i + 4)	r 2	-1 +	i + 2 2(i + 4)(i + 6)
	+ 2 Pe		Da m Da m + 4	2	i	w i	i 2 4(i + 2)(i + 4)(i + 6) + 2i + 8	r 2	-1 .	(B 2)
										4)
	+ 2 Pe	Da m Da m + 4		i	w i	1 i + 2	2	r i+2	-1
	-2	1 i + 4	2		r i+4	-1 +

† Email addresses for correspondence: yohan.davit@imft.fr, sylvie.lorthois@imft.fr

Acknowledgements

Part of this work has been done while S.L. was the Mary Upson Visiting Professor, Meinig School of Biomedical Engineering, Cornell University. The authors gratefully acknowledge P. Tsai, P. Blinder and D. Kleinfeld for sharing anatomical data, N. Nishimura and C. Schaffer for useful discussions and A. Ghigo for carefully reading the manuscript. They also gratefully acknowledge the anonymous referees for their constructive comments and suggestions. They thank the COSINUS department of IMFT for help on computational aspects. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 615102 (https://erc.europa.eu/). M.B. was the recipient of an international mobility grant from Ecole des Docteurs de Toulouse. This work was performed using HPC resources from CALMIP (grant no. 2016-P1541).

Supplementary movie

Supplementary movie is available at https://doi.org/10.1017/jfm.2019.866.

Appendix A. Effective transport coefficients

In § 2.3.2, we introduced the different classes of transport equations obtained using multiscale asymptotics in the regime of weak couplings. Each class of equation corresponds to a roman numeral in figure 4, and has a different structure. The transport equations associated with domains II, IV, V, VI, VII involve effective transport coefficients and we give their expression below.

A.1. Domains II and IV Equations (2.67) and (2.68) introduce an effective reaction rate, K eff which expression depends on the scaling of the membrane Damköhler number as follows:

where a is the closure variable defined by C V,1 = Da m a(ρ)C V,0 . This closure variable is the solution of the following equation ∇ 2 ρ a = -2, which verifies the boundary condition at the vessel wall n • ∇ ρ a(1) = -1.

A.2. Domain V Equation (2.59) introduces an effective diffusion coefficient, D eff . This domain corresponds to Taylor's dispersion regime and is treated in detail in § 2.3.1. Here, we have

• b where b is the closure variable defined by C V,1 = b(ρ) • ∇ z C V,0 and is the solution of the following equation Ũ = ∇ 2 ρ b, which verifies the boundary condition at the vessel wall n • ∇ ρ b(1) = 0.

A.3. Domain VI Equation (2.69) introduces an effective reaction rate, K eff , and an effective diffusion coefficient, D eff . For this domain we have

where b is the exact same closure variable as in domain V.

A.4. Domain VII Equation (2.70) introduces an effective reaction rate, K eff , an effective diffusion coefficient, D eff , and an effective velocity, U eff . For this domain we have