Reduced-order modeling of hemodynamics across macroscopic through mesoscopic circulation scales

TL;DR

This paper introduces a reduced-order hemodynamic model that bridges macroscopic and mesoscopic blood flow scales, accurately simulating vascular responses across different vessel types and scales.

Contribution

It presents a novel multi-scale modeling approach combining stochastic vascular geometries, scale-specific pruning, and diverse structural models for blood flow simulation.

Findings

Numerical results agree well with experimental data.

Model captures transition in flow pulsatility and wall shear stresses.

Effective in simulating microcirculation responses.

Abstract

We propose a hemodynamic reduced-order model bridging macroscopic and meso-scopic blood flow circulation scales from arteries to capillaries. In silico tree like vascular geometries, mathematically described by graphs, are synthetically generated by means of stochastic growth algorithms constrained by statistical morphological and topological principles. Scale-specific pruning gradation of the tree is then proposed in order to fit computational budget requirement. Different compliant structural models with respect to pressure loads are used depending on vessel walls thicknesses and structures, which vary considerably from macroscopic to mesoscopic circulation scales. Nonlinear rheological properties of blood are also included and microcirculation network responses are computed for different rheologies. Numerical results are in very good agreement with available experimental measurements. The computational model captures the dynamic transition between large-to small-scale flow pulsatility speeds and magnitudes and wall shear stresses, which have wide-ranging physiological influences.

Tables

Table 1. Table 1: Representative vessels distribution and size in the human circulatory system.[ 25 , 10 ] Here the number n 𝑛 n is typical of the chosen averaged diameter D 𝐷 D presented in the table. Vessels categories may as well be thought in terms of size ranges; for instance, small arteries in [ 600 − 100 ] ​ μ ​ m delimited-[] 600 100 𝜇 𝑚 [600-100]\mu m , arterioles in [ 100 − 10 ] ​ μ ​ m delimited-[] 100 10 𝜇 𝑚 [100-10]\mu m and capillaries in [ 10 − 5 ] ​ μ ​ m delimited-[] 10 5 𝜇 𝑚 [10-5]\mu m .

	macroscale			mesoscale			macroscale
	Aorta	Arteries		Arterioles	Capillaries	Venules	Veins
		large	small
Number $n$	1	40	7100	$1.6\cdot {10}^8$	$5.5\cdot {10}^9$	${10}^9$	${10}^2$
Lumen diameter $D$	$2.6cm$	$8mm$	$600\mu m$	$20\mu m$	$9\mu m$	$20\mu m$	$5mm$
Thickness $h(\%D$)	$5.8\%$	$8.2\%$	$18.6\%$	$32.5\%$	$10\%$	$\approx 10\%$	$\approx 10\%$

Table 2. Table 2: Numerical values used at each network bifurcation, during the stochastic generation of the vascular tree, to determine vessels diameter and length. Given a parent vessel of diameter D p subscript 𝐷 𝑝 D_{p} at a certain generation, daughter vessels diameters are obtained as: D d ​ 1 = κ ​ D p subscript 𝐷 𝑑 1 𝜅 subscript 𝐷 𝑝 D_{d1}=\kappa\,D_{p} , with κ ∼ 𝒰 [ κ min , κ max ] similar-to 𝜅 subscript 𝒰 subscript 𝜅 min subscript 𝜅 max \kappa\sim\mathcal{U}_{[\kappa_{\text{min}},\kappa_{\text{max}}]} and D d ​ 2 = ( 1 − κ m ) 1 m ​ D p subscript 𝐷 𝑑 2 superscript 1 superscript 𝜅 𝑚 1 𝑚 subscript 𝐷 𝑝 D_{d2}=(1-\kappa^{m})^{\frac{1}{m}}\,D_{p} , according to Murray’s power law. Concerning the branch length, except for the capillaries, the following deterministic relation with the diameter is used [ 26 ]: l . ​ ( D . ) = ζ ​ ( D . ) γ subscript 𝑙 . subscript 𝐷 . 𝜁 superscript subscript 𝐷 . 𝛾 l_{.}(D_{.})=\zeta(D_{.})^{\gamma} with m = d f + γ 𝑚 subscript 𝑑 f 𝛾 m=d_{\text{f}}+\gamma , and ζ 𝜁 \zeta determined according to the relations and measurements introduced in Kamiya et al. [ 26 ]. A fractal dimension of d f = 1.75 subscript 𝑑 f 1.75 d_{\text{f}}=1.75 is considered in this work.

$D_p$ range	${\kappa }_{\text{min}}(\%)$	${\kappa }_{\text{max}}(\%)$	$m$	$d_{\text{f}}$	$\zeta$	$\gamma$
$D_p\in [2.6cm,500\mu m]$	50	98	2.7	1.75	104	0.95
$D_p\in [500\mu m,100\mu m]$	80	95	2.7	1.75	1453	0.95
$D_p\in [100\mu m,10\mu m]$	65	79	3	1.75	1453	1.25
	75	95	3.5	n/a	n/a	n/a

Table 3. Table 3: Main physical and numerical simulation parameters.

Root vessel diameter $\left[mm\right]$	$D_{\text{root}}=1$
Minimal capillary diameter $\left[\mu m\right]$	$D_{\text{min}}=8.5$
Blood density $\left[kg\cdot m^{-3}\right]$:	$\rho =1050$
Plasma kinetic viscosity $\left[m^2\cdot s^{-1}\right]$:	$\eta ={10}^{-6}$
Discharge hematocrit:	$H_d=0.45$
Cross-section velocity profile:	$\alpha =4/3$ (parabolic)
Membrane Poisson ratio:	$\nu =1/2$
Membrane Young modulus $E$:	piecewise linearly interpolated from [43]
Venous pressure $\left[mmHg\right]$:	$p_{\text{venous}}=10$
Interstitial pressure $\left[mmHg\right]$:	$p_{\text{ext}}=10$
Mesh cells per vessel:	$e=1$
Legendre polynomial degree:	$p=7$
Over-integration factor:	$3/2$
Number Gauss-Legendre quadrature points/cell:	$n_q=11$
Time-step $\left[\mu s\right]:$	$\Delta t=0.3$
Simulation duration $\left[s\right]$:	$T_{\text{final}}=15$
Inlet BC (pulsatile half sinusoidal flow wave) $\left[ml\cdot mi{n}^{-1}\right]$:	${\overline{Q}}_{inlet}=1.74$
Systolic duration $\left[s\right]$:	$t_c^{\text{sys}}=0.3$
Pulse pressure frequency $\left[Hz\right]$:	$f=3.33$
Heart beat frequency $\left[bpm\right]$:	$f_{HB}=60$
Outlet BC (3-element Windkessel-type: RCR)	for $D_n\ge 100\mu m$
Outlet BC (1-element Windkessel-type: R)	for

Equations

$$\frac{\partial A}{\partial t}+\frac{\partial Au}{\partial x}$$

$$\frac{\partial u}{\partial t}+(2\alpha -1)u\frac{\partial u}{\partial x}+(\alpha -1)\frac{u^2}{A}\frac{\partial A}{\partial x}$$

$$p$$

$$k(\omega )=\pm \frac{\omega }{c_{\text{ref}}}{\left(\frac{1}{1-F_{10}(\omega )}\right)}^{1/2},\text{with}{F}_{10}(\omega )=\frac{2J_1(i^{3/2}{W}_{\alpha })}{i^{3/2}{W}_{\alpha }{J}_0(i^{3/2}{W}_{\alpha })},$$

$$D_p^m=D_{d1}^m+D_{d2}^m,$$

$$\eta =\left[1+({\eta }_{0.45}-1)\cdot \frac{(1-H_d{)}^C-1}{(1-0.45{)}^C-1}\cdot {\frac{D}{D-1.1}}^2\right]\cdot {\left(\frac{D}{D-1.1}\right)}^2,$$

$$C=\left(0.8+\exp (-0.0075D)\right)\cdot \left(-1+\frac{1}{1+D^{12}\cdot 10^{-11}}\right)+\frac{1}{1+D^{12}\cdot 10^{-11}},$$

$${\eta }_{0.45}=6\exp (-0.085D)+3.2-2.44\exp (-0.06D^{0.0645}),$$

$$f=-2\pi\underbrace{\mu_{\text{plasma}}\,\eta\big{(}D(A),H_{d}\big{)}}_{\mu}\frac{\alpha}{\alpha-1}u,$$

$$p=p_{\text{ref}}+\beta (\sqrt{A}-\sqrt{A_{\text{ref}}}),$$

$$p_{\text{ref}}$$

$$\beta$$

$$R_{\text{ref}}=\frac{R_{\Delta p_t=0}}{1-\frac{1-\nu }{E}{p}_{\text{ext}}},h_{\text{ref}}=\sqrt{R_{\text{ref}}^2+h_{\Delta p_t=0}(h_{\Delta p_t=0}+2R_{\Delta p_t=0})}-R_{\text{ref}}.$$

$$\frac{\partial \mathbf{U}}{\partial t}+\mathbb{H}(\mathbf{U})\frac{\partial \mathbf{U}}{\partial x}+\mathbf{S}(\mathbf{U})=0,$$

$$\mathbf{U}=\left(\begin{array}{c}A\\ u\end{array}\right),\mathbb{H}=\frac{\partial \mathbf{F}}{\partial \mathbf{U}}=\left(\begin{array}{cc}u& A\\ \frac{1}{\rho }\frac{\partial \psi }{\partial A}+(\alpha -1)\frac{u^2}{A}& (2\alpha -1)u\end{array}\right)\text{and}\mathbf{S}(\mathbf{U})=\left(\begin{array}{c}0\\ \frac{f}{\rho A}\end{array}\right).$$

$$\frac{d{\hat{\mathbf{U}}}_{i,{\Omega }_e^n}^j}{dt}=\mathcal{F}({\mathbf{U}}_{{\Omega }_e^n}^{\delta }),j=0,\dots ,p_e,i=1,2,e=1,\dots ,N_e,\text{and}n=1,\dots ,N,$$

$$\mathcal{F}({\mathbf{U}}_{{\Omega }_e^n}^{\delta })=-\frac{1}{J_{{\Omega }_e^n}}{\left(\frac{\partial {\mathbf{F}}_i({\mathbf{U}}_{{\Omega }_e^n}^{\delta })}{\partial x},L_j\right)}_{{\Omega }_e^n}-\frac{1}{J_{{\Omega }_e^n}}{\left[L_j\left[{\mathbf{F}}_i^u-{\mathbf{F}}_i({\mathbf{U}}_e^{\delta })\right]\right]}_{x_{{\Omega }_e^n}^{-}}^{x_{{\Omega }_e^n}^{+}}+{\left({\mathbf{S}}_i({\mathbf{U}}_{{\Omega }_e^n}^{\delta }),L_j\right)}_{{\Omega }_e^n},$$

$${\tau }_{rx}(x,t)\equiv \text{WSS}=\mu {\left(\frac{\partial \tilde{u}}{\partial r}\right)}_{r=R}=\mu (x,t)u(x,t){\left(\frac{\partial \varphi }{\partial r}\right)}_{r=R}\left(=4\mu \frac{u}{R}\text{ for Poiseuille profile}\right).$$

Full text

††articletype: Article Type

1 Introduction

Microcirculation plays a central role for many important physiological mechanisms of the circulatory system. It is the main site for exchanges between blood and tissues and also regulates and reorganizes the blood flow according to metabolic activity or the development of certain pathologies, e.g. neurodegenerative diseases, pulmonary arterial hypertension or tumors. Microvascular structure is highly complex but very well distributed across spatial scales, thanks to a fractal type of geometry organized in branching patterns. As mentioned by Secomb et al [1] “The microcirculation represents an important ‘mesoscale’ in physiological systems, functionally bridging higher and lower scales." Indeed the functional status of the microcirculation, which is associated from now on to the vascular mesoscale, is determined by complex processes taking place at the cellular/molecular micro- and nanoscale and thus strongly influences tissues and organs. Reciprocally, large macroscale systemic quantities such as blood pressure and flux balance inevitably affect the function of the microcirculation.

Despite tremendous progress, reliable in vivo measurements of microcirculation structure and function remain very difficult. Numerical models are needed to establish quantitative relationships with phenomena occurring on these larger and smaller scales. According to Secomb et al. [1] and Boyer et al. [2], several coupled and nonlinear submodels must interact at different scales in order to account for nonlinear blood rheology, microvascular morphology and extravascular pressure gradient. A blood flow model interacting with a microvascular architecture through vessels structural mechanical models and subject to a scale-dependent blood rheology model must be considered. Ideally, a solute transport model of blood compounds and its interaction with an arteriolar vasomotion model should also be considered. Nevertheless, for each of these components, it remains impossible to simulate phenomena across these spatial multi-scales by means of 3D models presented by Quarteroni et al. [3] One-dimensional reduced-order models (ROM) are commonly used to simulate blood flow regimes for which pulse waves propagate in large compliant arteries. They are increasingly being used to study physiological hemodynamic phenomena referenced by Mynard et al. [4] and to assist in the validation or measurement methods. Most studies focus on macroscales and investigate the impact of mechanical modeling and calibration of the arterial walls [5] , efficiency and robustness of deterministic and stochastic numerical schemes [6, 7, 8] and boundary conditions [9] or detailed anatomy [10].

While computational models of very different complexities and fidelities have been developed to model circulation at smaller scales,[11, 12, 13, 14] very few topologically/geometrically-detailed network-based attempts exist.[15, 16] Moreover, their reliability and use are still hampered mainly because they do not encapsulate all of the necessary submodels and are often dissociated from those of the systemic circulation. Other reasons include disregard of structural compliance and non-stationary flows, which are detrimental to the study of the significance of pulsatility in microcirculation revised by Pan et al. [17] and Lee et al. [18]. Indeed, several experimental studies have shown that the pressure pulse, generated by the heart, reaches the capillary or even beyond [19]. In fact, microvascular pulsatility has wide-ranging physiological influences, e.g. from the stimulation of nitric oxyde release [20] regulating the microvascular tone responsible for fine control of local blood flow to the control of inflammation response during cardiopulmonary bypass operations. Finally, non-linear blood rheology is often over-simplified in those models.

Computational domains representing large-scale complete and anatomically accurate microvasculature are very difficult to construct from biomedical images. Indeed, these highly distributed and space-filling structures bear a high degree of spatial heterogeneity in architecture and topology. The authors then rely on the inherited framework of one-dimensional ROM to simulate the blood flow, solved thanks to a discontinuous Galerkin (DG) method. In section (2), they show how to extend its capability toward microcirculation. For the computational domain, a method is proposed to generate idealized yet realistic vessel trees by means of stochastic algorithms. Several truncation strategies are also investigated in order to efficiently prune the tree and to provide calibrated boundary conditions in place of the missing branches. The vascular tree is mathematically described by a graph of blood vessels. Moreover, vessels arterial walls being very different in terms of their response to pressure loads from macroscopic to mesoscopic circulation scales, different structural models are adopted to address the specifics of the wall thickness-to-radius ratios. Finally, the authors have accounted for some of the blood rheology effects through the parametric calculation of a dynamic macroscopic description of the *apparent * blood flow viscosity. Some simulation results for a vasculature tree ranging from macro to mesoscale differing by more than two orders of magnitude are also presented and discussed in sections (3) and (4), respectively. It is through the inclusion of these several important multiscale modeling components, unified under a common computational framework capable of simulating the dynamic effects generated by a reasonably large circulatory network (i.e. $\mathcal{O}(10^{4})$ vessels) from large arteries to capillaries, that this work is original.

2 Materials and Methods

2.1 Microcirculation reduced-order model

One-dimensional ROM are commonly used to simulate convection-dominated blood flow for which pulse waves propagate in large elastic arteries.[21, 6, 7] They rely on a fluid-structure interaction mathematical framework that is much simplified when assumptions of Newtonian properties, linear elasticity (and to some extent nonlinear viscoelasticity) and homogeneous geometry are made for the blood fluid, the vessel walls mechanical response and the circulation network, respectively. They predict hemodynamic quantities with satisfactory accuracy and their low computational cost enables the simulation of broad arterial networks.[10] authors of the present paper rely on this inherited framework to enrich it and extend its capability toward microcirculation at mesoscopic scales both in terms of modeling and computational efficiency.

2.1.1 One-dimensional modeling framework

We consider an arterial tree of $N$ vessels, mathematically described by a graph: i.e. a network of thin, deformable, and axisymmetric arterial segments filled with blood, taken as an incompressible non-Newtonian fluid. The formulation directly derived from Navier-Stokes equations, for each straight arterial segment of lumen diameter $D$, based on the conservation of mass and momentum equations and a pressure law, takes the form of a nonlinear system:

[TABLE]

where $t$ denotes time, $x$ is the axial coordinate, $A(x,t)=\pi D(x,t)^{2}/4$ is the circular cross-sectional lumen area, $u(x,t)$ and $p(x,t)$ are cross-sectional averaged velocity and internal pressure, $\rho$ the blood density, $f$ the friction force related to the choice of axisymmetric velocity profile through $\displaystyle f=-2\mu\pi\frac{\alpha}{\alpha-1}u$, $\mu$ is the apparent blood dynamic viscosity, $\displaystyle\alpha=\int_{0}^{R}2r\,\tilde{u}^{2}/(uR)^{2}dr$ is the momentum flux correction (Coriolis) coefficient accounting for the nonlinear integration of radial velocities in each cross-section  [3], and $\tilde{u}=u\,\phi(r)$, with the velocity profile shape $\phi(r)=\frac{\alpha}{2-\alpha}\left(1-\left(\frac{r}{R}\right)^{\frac{2-\alpha}{\alpha-1}}\right)$.

The “tube law” relation, links the blood pressure to the vessel wall deformation, through a chosen function $\psi$, satisfying $(\frac{\partial\psi}{\partial A}>0,\;\psi(A_{\text{ref}})=0)$, which may reflects compliant effects through some elastic model, cf. subsection (2.1.4). The underscript $\cdot_{\text{ref}}$ denotes quantities in a reference unloaded state, $A_{\text{ref}}$ and $p_{\text{ref}}$ are taken both constant along each artery.

Due to the interaction between the pulsatile blood pressure and the vessels distensibility $\displaystyle\mathcal{D}=\frac{1}{A}\frac{dA}{dp}$, flow/pressure waves travel at a certain speed $c(x,t)$ in the system. Wave propagation and attenuation strongly depend on $\mathcal{D}$ and the ratio between fluid pulsation and viscous forces through the Womersley number: $W_{\alpha}=\frac{D}{2}\left(\frac{\rho\omega}{\mu}\right)^{1/2}$, proportional to the oscillation angular frequency $\omega$.

It is worthwhile reminding the main results obtained for a linearized version of Equation 1 with purely elastic arterial wall and viscous fluid under different regimes. The approach considers small perturbations and looks for solutions of the form: \big{(}u,A,p\big{)}(\omega,x,t)=\big{(}\hat{u},\hat{A},\hat{p}\big{)}(\omega,x=0)\exp{(i(\omega t-kx))} with $k=k_{r}+i\,k_{i}$, the complex wave number with $c=\omega/k_{r}$ the wave speed and $\gamma(\omega)=-2\pi k_{i}/k_{r}$ the attenuation coefficient.[22] Different scenario take place depending on the value of $W_{\alpha}$:

• –

for large $W_{\alpha}\rightarrow\infty$, the viscous term is neglected and simplification of the system leads to a pressure wave equation. The wave speed $c(w)=c_{\text{ref}}=\sqrt{1/\rho\mathcal{D}_{\text{ref}}}$ is constant with no wave attenuation: $\gamma(\omega)=0$. Concomitant forward/backward pressure waves travel at constant speed without damping.

• –

for small $W_{\alpha}\rightarrow 0$, inertial term is neglected and the viscous term is modeled by a Poiseuille law. This simplification leads to a pressure diffusion equation. The wave number is complex with: $c(w)=\frac{1}{2}W_{\alpha}c_{\text{ref}}$ and $\gamma(\omega)=2\pi$.

• –

for arbitrary $W_{\alpha}$, the wave number is:

[TABLE]

from which wave speed and attenuation are readily computable ($J_{0,1}$ are Bessel functions).

The computational model should be able to accurately capture all different pulse waves interacting at different speeds through the system. The construction of the computational domain, blood rheology model and vessel wall structural models will be described next. Let us start with the generation of a healthy vascular network.

2.1.2 Vascular tree generation and graph models

Microcirculation models of healthy tissues are usually established based on physical microvascular network structures,[17, 18] which, by nature, bear some degree of spatial heterogeneity in architecture (vessels diameters and lengths, vascular density) and topology (hierarchical organization, deviation from symmetry). Relative to normal tissues, tumoral tissues microcirculation (not treated here) exhibits higher structural and functional heterogeneity. With great effort, microvascular networks may be reconstructed from measurements (e.g. by intravital microscopy). The vasculature is organized in different patterns depending on the scale under consideration, e.g. arteries make a divergent arborescence, while the capillaries form some kind of “mesh-like" space-filling network.[23]

Another approach is to generate vessel trees ex novo by means of mathematical algorithms, across all scales down to the capillaries size, (statistically) constrained by some realistic morphological and topological principles, e.g. table (1).[24, 22]

Some notations and concepts

We now consider a circulatory tree $\mathcal{T}$ organized into small arteries, arterioles and capillaries (draining venules and veins are voluntarily not included and replaced by appropriate boundary conditions). The tree is split into $N$ vessels $v_{.}$, such that $\mathcal{T}=\bigcup_{n=1}^{N}v_{n}$ . Each vessel $v_{n}$ has its own constant length $l_{n}$ and reference lumen diameter $D_{n,\text{ref}}$ (we will often omit the $n$ subscript for convenience). Junctions between vessels are simply modeled as nodal points.

Biological branching systems most often exhibit fractal nature, i.e., a scale-independent self-similarity in the bifurcation pattern of their architecture. The measured fractal dimensions of vascular systems reveal a power function with a single exponent $<\cdot>^{-d_{\text{f}}}$ expression. This fractal dimension has been reported to be $d_{\text{f}}\in[1.1,1.86]$[23]. On the basis of a stochastic model relating the probability of branch density with aggregated branch length, morphological parameters such as surface area or content volume may be expressed through definite integrals of $d_{\text{f}}$ for a large range of diameters.[26] Here, diameter distribution is based on the law of Murray, which relates parent and daughter vessels at bifurcations as in Murray et al. [27]:

[TABLE]

with $m\approx 3$***Reported measured values in various vascular beds give $m\in[2.7,3]$ [24]. This universal law [28] was proposed for the optimal design of blood vessels, based on a balance between blood circulation and metabolic powers. While this law holds true for mid-scale circulation, this is not the case for the larger vessels, where inertia and turbulence effects may occur and for microcirculation due to shear-thinning effects. When Fåhreus-like phase separation effects occur in the latter case, thanks to some simplified blood constitutive model (e.g. Quémada’s model), it is possible to derive an extended version of Murray’s law, and see how it affects the optimal geometries of fractal trees to mimic an idealized arterial network. In this case, it was shown that the geometric scaling becomes less obvious as the optimal factors are now dependent on the generation index, the size of the root vessel and the relation between vessel diameter and mean shear rate [29]. In the following, a simpler hands-on approach is chosen, for which the Murray’s law is calibrated differently depending on the vascular scale.

For a given parent diameter $D_{p}$, a relation of proportionality $D_{d1}=\kappa\,D_{p}$ is assumed, and $D_{d2}$ is obtained from Equation 3. The branch length-diameter relationship follows a power law: $l_{n}=\zeta\times D_{n,\text{ref}}^{\gamma}$,[24] with $m=d_{\text{f}}+\gamma$, and $\zeta$ a scaling determined according to other relations, cf. Takahashi.[24]  

In our case, the coefficients $\kappa$ and $(m,\zeta,\gamma)$ will be random and deterministic, respectively. Their respective adjustments will depend on the vascular type and domain as explained next. As derived in Kamiya et al.[26], it is possible to relate the quantity $\zeta$ to the knowledge of the fractal dimension, the junction exponent, geometrical bounds of largest $R_{max}$ and smallest $R_{min}$ vessel radius in a fractal tree and its approximate volume $V$: $\zeta=\lambda\iota/R_{max}^{d_{\mathrm{f}}+\gamma}$ where $\lambda=V(-d_{\mathrm{f}}+2)/(\pi\iota(R_{max}^{-d_{\mathrm{f}}+2}-R_{min}^{-d_{\mathrm{f}}+2}))$ with $\iota=d_{f}/(R_{\mathrm{min}}^{-d_{\mathrm{f}}}-R_{max}^{-d_{\mathrm{f}}})$. The entire adult human arterial volume approximate $V\approx 700\mathrm{ml}$ from aorta down to capillaries [26 ]. Data provided by Boileau et al.[6] allow us to frame the volume of the macrocirculation tree. Under the assumption of a fractal tree, the corresponding $\zeta$ parameter can be found thanks to the model. The same procedure is then used to evaluate the microcirculation tree length parameter.

Graph representation and Stochastic growth model

For computational purpose, the topological network is reduced to a graph, described by node coordinates and a connectivity matrix. The graph is a full binary tree, which every node has either zero or two children. Its representation is very handy as it naturally provides (statistical) information related to the connectivity and the geometry of the entire (or part of the) network. Our tree generation is inspired by stochastic growth models (e.g. RSB model[30]). It starts from a root vessel ($D_{\text{root}}$) and adds new segments at each generation in an iterative stochastic process, until a stopping criteria is met (e.g. minimal capillary size). After numerous trials, our search of realistic trees, provided us with robust values of $(\kappa,m,\zeta,\gamma)$ as described in table (2).

The bounds of the $\kappa$ probability measure are adjusted in order to provide vessels distributions, whose statistics approach our reference, cf. table (1). We notice that the $\kappa$ range for the arterioles is consistent with the optimal prediction of Moreau & Mauroy [29]. We note that the Murray’s exponent has to be slightly adjusted for arteries and arterioles. While the same $d_{\mathrm{f}}$ is considered for the entire arterial tree, the calibration of $\zeta$ is made on two different scales, i.e.: on a fractal tree with extremum diameters $(12mm,500\mu m)$ and volume $200ml$,[26, 25] and on a more distal tree with diameters $(600\mu m,8\mu m)$ and volume of $0.043ml$.[26]

Concerning the capillaries we still generate them with a tree structure but, as the capillary bed is often made of a random homogeneous mesh-like structure, we have to increase the exponent to a larger value ($m=3.5$), in order to produce a larger number of vessels, in agreement with physiological studies [31]. Moreover, we rely on a random distribution to generate their lengths $l_{\text{cap}}\sim\mathcal{U}_{[500\mu m,700\mu m]}$ according to the literature.

Figure 1 shows an example of a generated network ranging from macroscopic ($D_{\text{root,ref}}=1mm$) to mesoscopic ($D_{\text{min,ref}}=7\mu m$) scale. The feeding vessel, despite being of moderate size (i.e. $1mm$), does represent the macroscale. Without restriction, it is possible to start with a larger vessel (e.g. $\mathcal{O}(cm)$), by connecting to the current model a network of large arteries such as the 55/56-artery networks frequently used in the literature, [32, 6, 7]. Not all vessels are represented in the figure, but the connectivity, length and diameter of the depicted segments are up to scale. Zoomed-up inset plots reveal some fine details (capillaries represented as red vessels). Figure 2 summarizes different geometric statistical markers along the different generations of that network, differentiating among small arteries ($D>100\mu m$, in green), arterioles ($10\mu m<D<100\mu m$, in blue) and capillaries ($D<10\mu m$, in red). The distributions in Figures 2(a)-2(b),2(d) exhibit very similar features as the ones described by imaging capabilities, such as: the spread of the capillaries across $\sim 20$ generations, a fraction of capillaries close to $60\%$, a hyperbolic decline of averaged arteriolar diameter with vessel generation, as in [33]). The topological heterogeneity of the network is also well captured in the statistics of Figure 2(d), with very different diameters (and lengths, not represented) within a given generation. Skewed distributions (with positive skewness) are obtained for small arterioles (here after the 17${}^{\text{th}}$ generation), in agreement with the literature [33]. These types of analysis were repeated several times for the generation of multiple networks. By doing so, it was checked that the proposed algorithm produced statistically robust and coherent graphs.

Truncation strategy

Hemodynamics simulations of such large-scale complex networks as in Figure  1 are obviously out of computational reach. We must propose a proper way of truncating the network and provide calibrated boundary conditions in place of the missing branches such that the simulation in the reduced domain produces a sound response. We have designed and tested several methods to prune the network based on deterministic or probabilistic trimming around one (or more) typical diameter(s) or generation scale(s) or based on imposing a decaying sequence of vessel growth inspired by the complete tree structure. The former approaches resulted in quite irregular distributions of vessels across generations. Moreover algorithmic probabilistic ingredients seemed essential in the selection process in order to avoid strong clustering. But these model reduction approaches remain useful if one is interested in keeping the full description of a particular region of the vasculature, while sacrificing large remaining domains through early-generation trimming.

The latter approaches gave best results. They rely on constructions that mimic the growth distribution of the complete graph. Figure 2(c) shows how the binary-type growth of the complete tree (green solid line) deviates from the exponential asymptote (oblique gray line) close to the seventeenth generation (cf. vertical green dashed line) due to the apparition of the capillaries, before to saturate and decline. This sequence is rescaled according to the targeted computational budget in order to keep a similar smooth profile (orange solid line) and probabilistically imposed to produce a more regular trimming. The new distribution departs earlier from the asymptote due to the imposed trimming. Nevertheless, the right tail of the distribution persists across a generation range over which capillaries take birth, insuring the presence of those (although in smaller number) in the truncated network.

Figure 3 shows an example of a truncated network generated according to that pathway and ranging from macroscopic ($D_{\text{root,ref}}=1mm$) to mesoscopic ($D_{\text{min,ref}}=8.5\mu m$) scale. There are 312 capillaries (not visible) quite uniformly distributed within the network.

We have now to include the blood rheological effects into the model.

2.1.3 Blood rheology

Low shear regimes exist in the blood circulation due to pulsatile, recirculating or bifurcating effects, making blood a shear-thinning, viscoelastic fluid with non-Newtonian rheology. Aggregation and disaggregation of red blood cells (RBC) in elongated «rouleaux» structures as well as their elastic properties and related dissipation mechanisms are mainly the cause of this complex behavior which has significant influence on various healthy and pathological physiological processes, particularly at microcirculation scale. Most accurate time-dependent non-Newtonian constitutive models include viscoelastic, thixotropic or shear-thinning effects (e.g. Maxwell model [34] or Oldroyd-B model [35]). Nevertheless, validation in large-scale arterial networks with vessel wall compliant effects are often neglected and only very few works are available for 1D applications [36].

In this work, we follow a simpler approach accounting for non-Newtonian effects through a viscous contribution. Indeed, as blood is made of a suspension of solid particles in a plasma fluid, it makes non-constant its apparent (or effective) viscosity. The volume fraction of solid blood components named hematocrit $H$ is almost entirely made of RBC. It is common to distinguish among instantaneous tube hematocrit $H_{t}$ and hematocrit flux passing though a point in a given time as discharge hematocrit flux, $H_{d}$. An interesting phenomenon appears in vessels of small size, e.g. $D<300\mu m$: a plasma layer (PL) forms at the lumen wall while RBC tend to concentrate closer to the vessel center. This heterogeneity is in fact at the source of complex nonlinear rheological blood properties in microcirculation. The literature highlights three main effects.[37, 38]

Given a certain $H_{d}$, numerous studies notice a decreasing trend in the apparent blood viscosity with decreasing vessel diameter ($D>20\mu m$) where it starts increasing significantly ($D\leq 20\mu m$). This causes the change in the apparent viscosity of a Poiseuille blood flow and is the so-called Fårhæus-Lindqvist effect.

The plasma skimming effect describes the difference between tube and discharge hematocrit, and therefore the asymmetric distribution of hematocrit in channels downstream a bifurcation. [39] Indeed, hematocrit transported near a junction, meet daughter branches of different sizes and flow rates. Considering the hematocrit and Poiseuille cross-sectional profiles, the distribution may become unequal. This effect is more significant for (feeding) vessels of diameter $<40\mu m$.

Finally, another effect results from the velocity difference between $H_{t}$ and PL within the vessels and how it affects the local blood flow. As $H_{t}$ velocity is much higher than PL velocity, the hematocrit profile varies along the vessel, inducing longitudinal variations of the PL speed. This is known as the Fårhæus effect.

Fårhæus-Lindqvist effect

In the literature, numerous in vitro/vivo works have collected and fitted their data to parametric empirical descriptions of the apparent blood viscosity. For in vivo data, they have proposed a phenomenological relation for the relative apparent viscosity $\eta$, ratio between the blood and plasma viscosities. Derivations details for these models, provided here, are given in these references [40, 41]:

[TABLE]

with:

[TABLE]

and:

[TABLE]

the nominal apparent viscosity for $H_{d}=45\%$. Therefore, in this work we account for blood rheology effects through a dynamic macroscopic description of an apparent viscosity, which impacts the friction term in the momentum equation balance 1:

[TABLE]

with $\mu_{\text{plasma}}$ taken as the water viscosity and $D$ in $\mu m$. We have implemented and successfully tested this model against some literature benchmarks [18]. We have verified that accounting for the Fårhæus-Lindqvist effect induced an increased systolic flow by 10 to 45% when reaching $30\mu m$ diameter arterioles, while decreasing it by 10 to 30% in capillaries, compared to Newtonian simulations. Those results supported the choice of considering such effect for $D<300\mu m$.

Fårhæus and plasma skimming effects

The use of a transport equation is a simple attempt to model (part of) the local Fårhæuseffect. In this case, $H$ is passively transported at the flow velocity and retro-acts on the system through the blood viscosity, cf. Equation 7. We have implemented and validated this approach within the DG framework (cf. 2.2.1) for stationary blood flow velocity in a compliant vessel. Nevertheless, as mentioned in, [42] for pulsatile or even reverse time-dependent flow, we have noticed some spurious compression effects arising in the hematocrit concentration. These oscillations were not only detrimental to the numerical simulation (unstable situation) but they also provided very different results compared to the much smoother observed concentrations for experimental injection of a passive compound in arterial circulatory subsystem. Due to these difficulties, we have decided not to include the transport equation in this study, to the price of disregarding the plasma skimming effect.

2.1.4 Vessel wall structural model

The vessel walls of arteries and arterioles are composed of distinct layers. For instance elastic fibers included in the tunica media of large arteries provide vessels with an elastic component. Dynamic changes in the vessels geometry reflect both active and passive regulation mechanisms. Transmural pressure $\Delta p_{t}$ passively affects the lumen diameter. Ignoring the shear and inertial components, this effect may be modeled thanks to a Laplace’s law, leading to an algebraic relation. This simple thin membrane law has been successfully exploited in many applications, but is not well adapted for arteriolar wall bearing a thick muscularis layer designed to sustain high transmural pressures. In the following, we adopt the structural models described in [43] in order to address the specifics of the wall thickness-to-radius ratios.

Each vessel segment is modeled as an elastic but incompressible cylindrical shell with its own rigidity and only moderate deformations are considered. Moreover, due to the lack of experimental data from the literature, Young modulus $E$ is considered constant for each vessel, through the range of transmural pressure. In practice, given a reference vessel diameter, $E$ is piecewise linearly interpolated from literature data.[43] It was checked that other types of interpolant (nonlinear polynomials, spline, shape-preserving piecewise cubic) did not significantly change the value of the Young modulus within the considered range (i.e. $\ell_{2}$-norm difference $\leq 7\%$).

In Figure 4, we report the notation required for the definition of the relevant configurations. We refer to $(R_{\text{ref}},h_{\text{ref}})$ to describe its reference undeformed condition. First are the experimentally measured in vivo geometries, then comes the unloaded configurations: stress-free geometry for thick-walled rings and, zero transmural pressure geometry for thin-walled rings, respectively. Finally, predicted geometries are described.

The vessel wall motion induced by the transmural pressure is described by a linear elasticity relation, based on the reference unloaded state:

[TABLE]

where the $\beta$ parameter may encompass a vector of parameters describing the mechanical properties of the membrane, related to the distensibility of the vessel: $\beta=2/\mathcal{D}\sqrt{A}$. Causin et al.[43] have identified two cases:

1. –

Thick-walled case: vessel wall thickness is significant, e.g. $h\sim 35\%D$, then:

[TABLE]

Geometries measured in vivo $(R_{\text{d}},h_{\text{d}})$ do not correspond to unloaded conditions. In practice, an inverse problem has to be solved, whose unknowns are $(R_{\text{ref}},h_{\text{ref}},p_{d})$. As we do not rely on medical imaging and measurements, we assume without loss of generality that the measured configuration corresponds to zero transmural pressure, i.e. $A_{\text{d}}=A_{\Delta p_{t}=0}$, $h_{\text{d}}=h_{\Delta p_{t}=0}$. Under this assumption, the problem may be inverted:

[TABLE] 2. –

Thin-walled case: i.e. $h<10\%D$, which often leads to the following implicit simplification: $h\ll D$. Expressions in Equation 9 simplify to: $p_{\text{ref}}=p_{\text{ext}}$ so the transmural pressure is null $\Delta p_{t}=0$, the reference surface becomes $A_{\text{ref}}=A_{\Delta p_{t}=0}$, $h_{\text{ref}}=h_{\Delta p_{t}=0}$ and $\beta=\frac{E\,h_{\Delta p_{t}=0}}{(1-\nu^{2})\sqrt{\pi}\,R_{\Delta p_{t}=0}^{2}}$.

In Figure 4(b), we show examples of transmural pressures obtained from the two different models. In the simulations, structural model will be selected based on vessel wall thickness, cf. table (1).

2.2 Numerical approximation method

Getting back to our original ROM describing a single vessel in Equation 1, we now write it under its quasi-linear form:

[TABLE]

with,

[TABLE]

Because $A\geq 0$, the Jacobian matrix $\mathbb{H}$ possesses two real and distinct eigenvalues of different signs: $\lambda_{1,2}=\alpha u\pm\sqrt{c^{2}+\alpha u^{2}(\alpha-1)}$, with c=\big{(}\frac{A}{\rho}\frac{\partial\psi}{\partial A}\big{)}^{1/2} the pulse wave speed, and Equation 11 is strictly hyperbolic. It is then possible to diagonalize the system and perform a characteristics analysis, which is very useful to treat the Riemann problems involved at vessels interface and boundary conditions.[3] In Section (2.1), we have briefly motivated the fact that the source term will have an effect of varying stiffness depending on the flow regime. Inspired by earlier approaches,[17, 18] authors of the present document have decided to capitalize on the numerical solver documented in our previous studies and inherited for the case of convection-dominated macroscale arteries. [5, 7] Instead of domain-decomposing the network, in order to select most efficient numerical schemes depending on flow regimes, we have instead treated the network as a whole and computationally improved our solver to handle the various blood flow scales of our system. As explained later, much efforts have been dedicated to time-stepping tuning and memory constraints originating from the microcirculation.

2.2.1 Discontinuous Galerkin scheme

We adopt a DG method with a spectral$/hp$ spatial discretization with $N_{e}$ non-overlapping regions in each vessel and Legendre $L_{j=0,\dots,p}(x)$ polynomials, which is a very efficient framework for high-order discretization of convection-dominated flows, with low dispersion and diffusion errors. The approximated solution is noted: $\mathbf{U}^{\delta}$ in the network domain $\Omega$, and $\mathbf{U}^{\delta}_{\Omega_{e}^{n}}$ refers to its components in a particular cell $(e)$ of one of its arterial segment $(n)$. The discretized problem to solve forms of a system of ordinary differential equations for the modal coefficients of $\mathbf{U}$:

[TABLE]

where the right hand-side is made of several contributions:

[TABLE]

with inner products evaluated by numerical integration, thanks to Legendre-Gauss-Lobatto quadratures. An accurate calculation of the upwinded fluxes $\mathbf{F}^{u}_{i}$ across inter-elemental boundaries is required to solve the Riemann problem arising at each element interfaces. This is achieved by an efficient Riemann solver of Roe.[21] This well documented scheme has proven its efficiency for simulating hemodynamics in large- and medium-sized pulsatile arteries.[3]

2.2.2 Boundary conditions

The system requires a single inflow boundary condition (BC), and numerous outflow terminal BCs which are critical to the blood perfusion. A canonical time-dependent total volume flow rate is prescribed at the entrance of the network root in the form of: a half-sine function spanned over the systolic time (1/3 of the cardiac cycle) and null over diastolic time. At the exit of peripheral vessels, linear lumped parameter (0D) models are coupled to the global ROM, cf. for details about standard lumped parameter outflow models [32]. Matched three-element Windkessel models are used for small arteries and arterioles of size $D\geq 100\mu m$ while single-resistance models are used for smaller vessels and capillaries where fluid resistance dominates over wall compliance and fluid inertia.[16] While other BC closures are possible,[9, 44] we have estimated 0D BC parameters, i.e. resistance and compliance (when needed), from our network scaling laws (cf. table 2), under assumption of Poiseuille law and taking into account the effective non-Newtonian viscosity. This approach resulted in peripheral resistance and compliance distributions in agreement with the literature. For instance, we have recovered the exponential growth of the peripheral resistance with the increasing number of generations of bifurcations.[25] Single-resistance models attached to capillaries exit are tuned to obtain an average pressure drop of about $5\,mmHg$ over the tube.[14] In the three-vessel splitting flow junction, encountered many times in the system, the type of interface condition presented by Quarteroni et al.[3] which is based on consideration of the energy inequality in the vascular system, was chosen. We therefore impose mass and total pressure conservation. We mention that the current junction model does not account for momentum loss due to red cells and vessels deformation, bifurcation angles, mechanical properties, etc… but this approximation is probably less important in the microcirculation regime because of the dominance of the viscous forces and the resulting slower flow velocities.

3 Numerical simulations on an idealized human multi-scale vascular network

The new implementation of the numerical scheme allows dynamical simulations of large compliant networks. We have simulated fifteen cardiac cycles in a compliant vasculature ranging from macro to mesoscale differing by more than two orders of magnitude, cf. Figure 3. The choice of the different physical and numerical parameters are summarized in table (3). At the inlet, we have imposed an idealized half-sine input flow wave of period $t_{c}^{\text{sys}}=0.3s$, corresponding to the systolic ejection, for a total cardiac cycle of 1s duration. Two simulations were compared: one with constant blood viscosity and another one accounting for non-Newtonian effects through the dynamic evaluation of the effective viscosity, one of the objectives being to quantify some differences between the two time-dependent simulations.

3.1 Blood flow distribution

The entire tree was well perfused according to the mean volumetric blood fluxes observation, cf. Figure 5. The simulation generated an unsteady flow rate with an average of $1.74\,ml\cdot min^{-1}$ at the main feeding artery and $5\,nl\cdot min^{-1}$ in capillaries. The agreement between the data fit and experimental measurements in humans is excellent. [45] The decline of the corresponding averaged blood flow velocities (not represented) at three different scales along the arteriolar tree are also in good agreement with the literature, e.g. [22, 14]: for – small arteries at $D\sim 600\,\mu m$: $u_{\text{sa}}\approx 25\,mm\cdot s^{-1}$, – arterioles at $D\sim 40\,\mu m$: $u_{\text{art}}\approx 3\,mm\cdot s^{-1}$ and – capillaries at $D\sim 10\,\mu m$: $u_{\text{cap}}\approx 0.85-1\,mm\cdot s^{-1}$, respectively.

3.2 Pressure and pulsatility

The boxplot statistics of averaged (over space-time) pressure signals, plotted against vessel diameters, cf. Figure 6(a), show a strong pressure drop from the main feeding artery to the capillaries. As expected, mean pressure variability is largest across the arteriolar span, i.e. $100\mu m\leq D\leq 10\mu m$. Newtonian results (not presented) exhibit a stiffer pressure drop with significant prediction different far upstream in the network, as observed in Perdikaris et al.[44]. The dynamic characteristics of the pressure signals, captured by the pressure pulsatility index (here measured at the center of each vessel): $p_{I}=(p_{\text{systolic}}-p_{\text{diastolic}})/p_{\text{mean}}$, cf. Figure 6(b) also present strong decay with one order of magnitude difference from central to most distal regions. Nonetheless, pulsatility index in capillaries region do remain significant. Interestingly, statistical linear regressions of the first two $p_{I}$ moments against diameters (not included here) reveal a close-to-linear behavior of its mean ($r^{2}=0.98$; $p{\text{-value}}=1.8\cdot 10^{-9}$) and its standard deviation distribution ($r^{2}=0.96$; $p{\text{-value}}=7.8\cdot 10^{-7}$), respectively.

3.3 Wall shear stress and pulse wave speed

While our model contains simple non-Newtonian effects compared to other more complex constitutive models [36], it does include temporal and spatial structure-dependent viscous effects. In low shear regions, it is therefore important to quantify these effects. We monitor the wall shear stress (WSS), cf. Figure 7(a), that we relate to the wall shear rate through the effective viscosity $\mu(x,t)$:

[TABLE]

Only mean values over space/time are reported. Blue dots refer to the WSS obtained from the simulation with non-Newtonian rheology. As expected, despite a wide spread, WSS increase hyperbolically as arteriolar diameter decreases. Moreover, it is observed that WSS values are in good agreement with in vivo measures from recent literature.[47] Grey dots depict the Newtonian WSS response. They present values obtained from the simulation with constant blood viscosity ($\mu_{\text{constant}}=3.4\,Pa\cdot s$). The behavior is completely different: WSS increase stagnates for diameters lower than $50\mu m$. Non-Newtonian effects are therefore significant in most of the domain (not in small arteries) and must imperatively be accounted for.

The speed distribution of the pulse waves in the network are also worth investigation. Microvascular pulse wave velocity (PWV) emerges as a novel indicator of hypertension, but its dependence on various vascular scaling properties remains unclear. The theoretical PWV is evaluated , for a pressure pulse of frequency $f=3.33Hz$, thanks to Equation 2, under the favorable assumption whereby an antegrade pulse is transmitted free of reflection at branch points, and compared to the mean simulated wave speeds, cf. Figure 7(b) where only internal network vessels results are depicted. PWVs are directly processed from the time-dependent pressure and velocity wave signals. The PWV range is very wide between fast entering waves ($\sim 4m/s$) to slow travelling waves ($\mathcal{O}(cm/s)$) in the capillaries. The results agree well with the available literature. For instance, values obtained from the pressure signals in arterioles $D\approx 200\mu m$ are within a $[11-13cm/s]$ range compared to $[7-11cm/s]$ obtained by Intaglietta et al. [48] for the cat omemtum at slightly lower frequency $f=2.5Hz$, phase velocities increasing with frequency [19].

4 Discussion

The specifics of proposed modeling are discussed in details. The model bridges a macroscale $\mathcal{O}(cm-mm)$ to mesoscale $\mathcal{O}(\mu m)$ vascular territory and may be used to explore its baseline structure and function relationship through numerical blood flow simulations. Only arteries, arterioles and capillaries are modeled, avoiding the need for complex venular tree modeling.[43] Despite spatial dimension simplifications, hemodynamics simulation across several scales is a very demanding computational task. The new implementation makes feasible the dynamical simulation of large compliant network with $\mathcal{O}(10^{4})$ vessels, without resorting to domain decomposition, thanks to the rapid convergence properties exhibited by the spectral elements method, with reasonable dimensionality: spatial number of degrees of freedom: $2\times N\times e\times\lceil\frac{3p}{2}\rceil\sim\mathcal{O}(10^{5})$. Computer time to complete one cardiac cycle (with $\Delta t=0.3\mu s$) was 5 hours for an Intel Xeon E5 CPU having 3.7GHz clock frequency and 16Gb DDR3 of RAM. In case of larger systems, DG methods compactness (elements are discontinuous and inter-element communication is minimal) make them great candidates for efficient parallelization.[44] However, (semi-)implicit formulations are preferred for 3D models, the solver relies on an ($2^{\text{nd}}$-order Adams-Bashforth) explicit scheme. It performs well under convective stability constraint of CFL-type , i.e. $\Delta t^{\text{CFL}}\propto\frac{1}{2}\frac{\Delta x}{\lambda_{\text{max}}}$ ($\lambda_{\text{max}}$: maximum eigenvalue of the Jacobian matrix), but for system with large physical or numerical diffusive effects (low Reynolds and Womersley numbers), it is more strongly impacted by the diffusive stability constraint $\Delta t^{\text{diff}}\propto\frac{\Delta R^{2}}{\nu}$. Therefore, depending on the discretization retained, time-stepping choice insuring stability is not a straightforward task because of the interplay between these different effects. In Figure 8, we present an example of optimal time-stepping which is numerically determined for different discretization sizes. When axial discretization is fine enough, the time-step must asymptotically evolve as $\Delta t^{\text{CFL}}$, while at the opposite, it is fully dictated by the $\Delta t^{\text{diff}}$ scale, which is particularly stringent for vessel of smaller diameter.

In order to make the time-advancing computationally tractable, we have vectorized the computation of Equation 13, in the form of a dense matrix structure encapsulating the entire network information. Other studies have proposed to make implicit the time-integration, to the price of additional numerical effort in terms of improvement of system conditioning or sparse algebra acceleration [18, 17]. Another option, in the context of partitioned schemes, implemented in sequential or parallel fashion, is to resort to multi-level time stepping techniques,[49] where the different geometric components are solved with individual convenient time steps and matching interface conditions are fulfilled thanks to suitable synchronization procedures.

The model is conveniently compatible with a mixed bag of BCs and wall structural models, allowing fine tuning of physiological effects depending on the vessel scale. Moreover, the large scale representation authorizes a certain resilience to the choice of inaccurate outflow boundary conditions. Indeed, for macroscale systemic circulation, this is usually a difficult and unreliable task based on the calibration of lumped parameters models, that is somewhat mitigated in our case to the price of specifying the mechanical properties of the mesoscale vessels.

It is time-dependent and therefore offers great potential to study biomechanical determinants of pulsatility damping in great details through the vasculature. Results from Figure 7(b) quantify to what amount PWVs are slowed down along the arteriolar region. Accordingly with linear pulse transmission literature, pressure pulsation persists throughout the microvasculature and pulse propagation speed in small vessels is observed to be orders of magnitude less than predicted by the Moens-Korteweg relation. Despite the large variation, especially at large scale, potentially related to reflections at branching junctions, the agreement with the theoretical prediction is quite satisfactory [19]. Interestingly, pressure (velocity) wave speeds are a bit slower (faster) than the reference curve for medium size vessels, respectively. Eventually both quantities reach very comparable PWVs once they get closer to the capillaries, which is coherent with the assumption of synchronization of pressure and flux waves in microcirculation. More realistic viscoelastic model describing vascular wall properties would be interesting to include in our multiscale model, but should not qualitatively affect the results, despite increasing pulse attenuation and affecting pulse transmission time.

The model is applicable to large-scale vascular networks. A noticeable effort has been made in order to craft vascular structure with realistic topology. Some adjustments to the dyadic structure have been done in order to account for the capillaries density, and scale-specific pruning gradation of the tree has been proposed to minimize the sensitivity of the cut-off radius of the tree. . Nevertheless the effect of the tortuous morphology of the small vessels is neglected by the model. The arrangement of capillaries also strongly depends on the host organ and type of tissue. For instance, a network like configuration of the capillary bed, more realistic to model tumoral tissues, would require adding network anastomoses to the geometry (not a blocking point). Moreover, tumoral vasculature is made of tortuous and dilated vessels which are spatially strongly heterogeneous. In this case, Murray’s law does not hold anymore. But the proposed computational approach is completely general with respect to more complex and accurate large-scale human morphology, obtained from micro-computed tomography imaging for instance.[15]

In order to gain reliability and versatility, ROM must propose modeling improvements leading to significant error prediction reduction. Otherwise, low gains are very difficult to detect below the inherent approximation errors, and demanding solver developments become in this case hard to motivate. An experimentally-characterized Fårhæus-Lindqvist effect, inducing a diameter temporal dependence of viscosity, was added to the model and had a major impact on the overall flow, especially at the mesoscale. In order to account for phase separation effects, more involved modeling including an additional conservation equation for the RBCs may be introduced relating $H_{d}$ and $H_{t}$ through a time-dependent passive scalar transport equation where $H_{d}$ is advected at the flow velocity. We have implemented this approach but obtained unphysical oscillating hematocrit concentration results. These oscillations were not only detrimental to the numerical simulation (unstable situation) but they also provided very different results compared to the much smoother observed concentrations for experimental injection of a passive compound in arterial circulatory subsystem.[50]

While it is known that high-order DG methods without slope limiters dealing with problems governed by hyperbolic conservation laws suffer from spurious oscillations in the vicinity of discontinuities, we believe that more work is needed to understand the limitation of the one-dimensional cross-section averaged modeling for multi-scale transport. Constant concentration cross-section profiles are in this case unrealistic. Thus, uneven concentration of agent propagation in the vessel lumen and diffusion effects should be integrated in the ROM.

In conclusion, this ROM model may serve as a useful tool for various purposes, such as the investigation of pulsatility dependent microvascular physiology. Initial results have been encouraging and dynamical effects related to the coupling between the non-Newtonian blood rheology and flow pulsatility have been well captured. It may also serve as a foundation for the investigation of many pathologies related to microcirculation, which is crucial in the development of tumor-induced cancer but also in its detection and treatment, since it is the main path for the delivery of imaging agents and therapy drugs. Future works will involve model improvements for a better understanding of the micro-environment of vascularized tumors in relation to the systemic circulation. In particular, precise modeling of microvascularization is a major challenge today to strengthen and guide the interpretation of perfusion imaging used to monitor therapeutic effects.

Conflict of Interest

The authors declare that there is no conflict of interest regarding the publication of this article.

Acknowledgment

This work has been supported by the LabeX LaSIPS (ANR-10-LABX-32) managed by the French National Research Agency under the “Investissements d’avenir" program (ANR-11-IDEX- 0003-02).