{'fraction_non_alphanumeric': 0.03328284947904312, 'fraction_numerical': 0.01065635092969363, 'mean_word_length': 4.727244226146933, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mathematical modeling of local perfusion in large distensible microvascular networks

Microvessels -blood vessels with diameter less than 200 microns- form large, intricate networks organized into arterioles, capillaries and venules. In these networks, the distribution of flow and pressure drop is a highly interlaced function of single vessel resistances and mutual vessel interactions. In this paper we propose a mathematical and computational model to study the behavior of microcirculatory networks subjected to different conditions. The network geometry is composed of a graph of connected straight cylinders, each one representing a vessel. The blood flow and pressure drop across the single vessel, further split into smaller elements, are related through a generalized Ohm's law featuring a conductivity parameter, function of the vessel cross section area and geometry, which undergo deformations under pressure loads. The membrane theory is used to describe the deformation of vessel lumina, tailored to the structure of thick-walled arterioles and thin-walled venules. In addition, since venules can possibly experience negative transmural pressures, a buckling model is also included to represent vessel collapse. The complete model including arterioles, capillaries and venules represents a nonlinear system of PDEs, which is approached numerically by finite element discretization and linearization techniques. We use the model to simulate flow in the microcirculation of the human eye retina, a terminal system with a single inlet and outlet. After a phase of validation against experimental measurements, we simulate the network response to different interstitial pressure values. Such a study is carried out both for global and localized variations of the interstitial pressure. In both cases, significant redistributions of the blood flow in the network arise, highlighting the importance of considering the single vessel behavior along with its position and connectivity in the network.


Introduction
Microscopic blood vessels play the vital role of locally perfusing single body's organs. These circulatory districts include thousands of microvessels (diameter less than 200 µm), categorized as arterioles, capillaries or venules, according to their structure and function. First elementary assessments of microcirculatory mechanisms on skin or superficial organs date back at least to the 18 th century [1]. At present, techniques like positron emission tomography, magnetic resonance imaging and contrast echography allow to study in a non-invasive manner regional blood flow in internal organs of human patients. The information content of such measurements is, however, far to be complete, both in baseline conditions as well as in altered conditions. As a matter of fact, data such as vessel geometry, fluid-dynamics and physical parameters in physiological conditions or in altered conditions/provocation studies often cannot be coherently and comprehensively collected. These facts prevent from observing in detail the distribution of blood flow and pressure drop within the large and intricate microvascular networks, which are highly interlaced function of vessel resistances, in turn determined by blood biophysical properties, vessel mechanical and geometrical features [2,3]. For this reason, theoretical and computational models can help reducing this gap in knowledge.
Mathematical models in hemodynamics are present in literature since the 1960s. Nowadays, they are well assessed tools in the simulation of a limited number of vessels of major size like the aorta and collaterals, possibly coupled with reduced-order models for the rest of the circulatory system (see, e.g., the review works [4,5]). In the context of microcirculation, the presence of an exceedingly large number (> 10 4 ) of connected vessels with complex behavior makes the problem pretty different and requires to adopt specific techniques. Typically, models of microcirculation use sets of representative segments to describe a network of vessels of different size (for example regrouping large/small arterioles and venules, see, e.g., [6,7,8,9]). This approach maintains a low number of unknowns and allows to explore in a simple manner different regulatory mechanisms via phenomenological relations. What is lost is the spatial distribution of field variables, so that relevant geometrical and physical heterogeneities of the network cannot be represented, as well as complex internal interactions (see [10] for a discussion on this topic). Spatial heterogeneity is taken into account in a different set of papers, which represent relatively small vessel graphs as a collection of 1D distensible tubes, derived from mathematical algorithms [11,12,13,14] or from medical imaging data [15]. While these models include the effect of the geometrical localization of each vessel in the network the mechanical description is absent (rigid vessels) or very simplified. In this latter case, phenomenological vessel compliance laws are often used, which reproduce selected structural behaviors. Moreover, it is not kept into consideration the fact that microcirculatory networks, characterized by low values of the luminal pressure, comparable to the surrounding interstitial pressure, can experience severe reductions of the luminal cross section, till collapse, both in physiological conditions and -more dramatically-in presence of pathologies. The very few works addressing this phenomenon locate themselves at two opposites. Either, they study a single vessel (or very small networks) with a complex description, possibly 3D and anatomically accurate, [16,17,18,19,20,21], or they model mid-sized/large microvascular networks including a Starling/collapsible components [22,23,7] or phenomenological variations of the physical parameters (see [24] in a slightly different context). While the first approach is unable to be applied to a large network due to its huge computational cost, the second one is much more efficient but it does not consider the sophisticated coupling between vessel mechanics and pressure loads. These facts limit the spectrum of phenomena which can be analyzed.
In this article, we propose a mathematical model of general microcirculatory districts which is capable of dealing with large, general networks with an affordable time of resolution. Upon extracting the geometry from medical imaging data or from a computer-generated structure, each network is described as a graph of distensible tubes (vessels), representative of realistic structures under physiological pressures. A simplified fluid-structure interaction approach is considered. Namely, blood flow in each vessel is accounted for by a generalized Poiseuille's law, featuring a conductivity parameter which is function, among the others, of the area and shape of the tube cross section. To dispose of these latter parameters, we adopt thick or thin-wall structural models according to the physiological vessel wall thickness-to-radius ratio [25,26]. In addition, a buckling model for thin-walled vessels derived from [27] is included. Buckled vessels lose the circular-shaped cross section and assume an elliptical or dumb-bell configuration, causing a strong increase of vessel resistance to flow. The model can withstand partial or total vessel blockage, providing detailed information about the fluid-dynamics. This contrasts with simpler Starling resistor elements, which due to their switch-like behavior, cannot represent intermediate regimes (partial patency to flow). Using the network geometries proposed in [28,29], we simulate blood flow in the retinal circulation, which we already studied without including vessel compliance in [30,31]. The network is composed of more than 9000 vessels, with a tunable degree of asymmetry. After validating of the model against experimental measures (data from [32]), we carry out two sets of studies: i) we globally increase the external pressure, reaching the conditions for buckling to occur. We observe that the luminal pressure gradually increase along all the network, till buckling, after which a discontinuity in the behavior takes place, with a much more marked pressure increase and flow redistribution; ii) we locally increase the external pressure inside a spherical region located in correspondence of a region of the post-capillary veins. This perturbation is observed to extend its effect till four or five vessels generations away, with important redistribution of flow and resistance. An interesting, noticeable, key point emerging from the above results is the importance of vessel localization in the network. Vessels with the same mechanical and geometrical properties but laying in two different regions of the network display a pretty different behavior due to their local pressure levels and interaction with other vessels.
The paper is organized as follows: in Sect. 2 we present our mathematical model for microcirculatory districts. Namely, in Subsect. 2.2 we introduce the mathematical model for blood flow in a single vessel and the generalized Ohm's law connecting flow rate and pressure drop via the conductivity parameter. This latter is obtained in a coupled manner from the wall structure models as discussed in Subsect. 2.3 in pre-buckling conditions (see 2.3.1), where a particular attention is devoted to the Young modulus choice, and buckled conditions (see 2.3.2). In Sect. 2.4, we report a summary of the computation of the conductivity parameter. In Sect. 2.5, we discuss the importance of using a correct unloaded configuration, describing the numerical technique applied to compute it from measurements. In Sect. 3, we introduce the nomenclature to deal with a network and we present the conditions to couple single vessels converging in a node. In Sect. 4, we provide a summary of the model (see 4.1) and we discuss the numerical procedure employed to discretized the fully coupled problem (see 4.2). Then, in Sect. 5, we first introduce the network geometries we will use in simulations (see 5.1) along with the physical parameters (see 5.2) then, we present the results of the numerical simulations in different test cases (see 5.3). Eventually, in Sect. 6, we draw the conclusions, discussing the results along with their significance, the limitations of the model and the forthcoming work.

Microcirculation model 2.1 Geometrical model
The geometry, denoted in the following as the "measured geometry", of the network can be originally derived from digitized images or can be constructed in silico on a computer on the basis of anatomical data. In any case, our starting point is the 1D skeleton of the network along with its topological connectivity and cross sections and lengths distribution. Each segment of the skeleton represents a blood vessel and can branch at nodal junctions. To increase the computational accuracy, we introduce further subdivisions into elements on each segment (see [18] for a similar approach). We establish on each element a local system of cylindrical coordinates and we let the z-axis coincide with the element axis, arbitrarily choosing its orientation. The vessel element is endowed of the 3D structure of a straight cylinder of axis z, with uniform, but not necessarily circular, cross section. The r and θ coordinates lay in the plane of the vessel cross section (see Fig. 1). Elements belonging to the same vessel share homogeneous mechanical properties. From this geometrical model, we proceed by computing a reference ("unloaded") configuration as described in Sect. 2.5. This latter geometry represents the mathematical domain of the present model.

Blood flow model
The domain occupied by blood inside the vessel element (luminal space) is defined as (see Fig. 1) where A = (0, R(θ)) × (0, 2π), R(θ) being the position of the blood-wall interface and where L e is the element length. Blood circulation is modeled as the steady unidirectional flow of a Newtonian incompressible fluid with dynamic viscosity µ. Figure 1: Schematic representation of a portion of a microcirculatory network. Each vessel is described as a duct with straight longitudinal axis and is further partitioned into a series of consecutive short elements of arbitrary but constant cross section shape along the axis. The central part of the figure depicts one of such elements, of constant length Le, with highlighted the domain Ω f occupied by the blood flowing inside the luminal space and the wall structure. The right part of the figure represents the cross section A of the luminal space along with the thickness h of the vessel wall. The wall is internally loaded with pressure p l from fluid actions and externally loaded with given interstitial pressure pe. The local system of coordinates (r, θ, z) on the element is represented as well.
Letting p be the fluid pressure and u the axial (and only) velocity component, the continuity and momentum balance equations read: find p and u such that where ∆ rθ (·) is the Laplacian operator with respect to the (r, θ) coordinates. No-slip conditions are considered on ∂A × (0, L e ). Notice that the pressure field resulting from Eqs.
(1) has a constant gradient in the axial direction. Moreover, the pressure is uniform on each section, so that, straightforwardly, the fluid pressure p l acting on the internal surface of the wall structure is equal to p.

Remark 1.
In this work, we consider a dynamic viscosity depending, among the others, on the vessel cross section diameter (more generally speaking on the hydraulic diameter), as discussed in detail in Sect. 5.2, yielding a nonlinear coupling with the geometry. This aspect is dealt with numerically by resorting to an iterative technique, which results at each iteration the viscosity to be a constant, given, value, computed from quantities known from the past internal iteration (denoted here tout-court by µ with a slight abuse of notation).
Our goal is to obtain a form of Eqs. (1) which is amenable to be efficiently coupled with wall structure equations in the context of a large network of vessels. Introducing the non-dimensional variables r * = r/ R and u * = (−µ R 2 dp/dz −1 )u, R being a characteristic linear dimension of the cross section, we write the dimensionless form of Eq. (1) 2 as [33]: where A * = A/ R 2 , and u * = 0 on ∂A * . We use the solution u * of (2) to define the conductivity parameter [27] so that the expression of the volumetric flux of fluid can be re-formulated as the generalized Ohm's law connecting flux and pressure gradient Q = −σ dp dz .
Notice that the conductivity σ is a function of the geometry of the vessel cross section, this latter being itself an unknown of the problem. The coupling with a structural model for the vessel wall through the pressure loads closes the problem.
We now go back to the original fluid balance equations (1) and we integrate Eq. (1) 1 on the cross section area. Gathering the resulting equation and the Ohm's law (5), we obtain the (equivalent) system: find Q and p such that Remark 2. In system (6), the conductivity parameter must be constant in each vessel element. For this reason, albeit the blood pressure in each vessel element is a linear function of the z coordinate, we consider the structure to be loaded on the lumen interface with a unique constant pressure p function of p (for example, its average along the element length) and, thus σ = σ(p). This approximation is acceptable if the number of elements in each vessels are chosen in a such a way that the pressure gradients are not excessively high.

Vessel wall model
In order to compute the vessel conductivity from Eq. (3), we need to dispose of the vessel cross section area and shape as a function of the pressure loads. In other words, we must build via a structural model a tube law, mathematically connecting the vessel cross section area with the transmural pressure p t , defined as the difference between the luminal pressure p and the interstitial pressure p e [34]. We anticipate in Fig. 2 the tube law resulting from the present model. Observe, in particular, the different behavior of arterioles (thick-walled vessels) and venules (thin-walled vessels). Observe also how, for these latter, there exists a physiologically plausible value of transmural pressure under which the tube is not any more circular but assumes a buckled configuration. Notice that the cross section does not need to be completely closed for the vessel to be "functionally lost to the network". As a matter of fact, it suffices the section to be small enough to prevent red blood cells passage to compromise its physiological function [15]. In the structural model, each vessel element is modeled as an elastic ring made of elastic (Young modulus E) and incompressible material (ν=0.5), assumed to be circular in undeformed conditions (radius R u , thickness h u ). The same cylindrical coordinate system of the fluid model is considered, if not otherwise specified. Small deformations are considered, similarly to several works in this field, see e.g., [35,4]. In Fig. 3, we report the notation required for the mathematical discussion and the definition of the relevant configurations we will consider.

Structural model for pre-buckling transmural pressure
On applying the internal and external pressure loads, radial and circumferential stresses arise in the ring. We assume axisymmetry and plane stress conditions. Let η = η(r) be the radial displacement of a point of the vessel wall. Then, we establish the strain-displacement relations (refer also to Fig. 4) ε N being the radial strain and ε T the circumferential strain, respectively, and the pseudoelastic constitutive equations σ N being the principal radial stress and σ T the principal hoop stress, respectively, and E = E(p t ) a functional representation of the Young modulus to be discussed at the end of this section. We close the problem considering the equilibrium equation Figure 3: Characteristic configurations in the structural models. Configuration I is the experimentally measured "in-vivo" geometry, supposed to be circular. The arrows indicate the steps followed in the computations to obtain a certain configuration from this configuration. According to the different modeling chosen as a function of the wall thickness-to-radius ratios (see Sect. 2.3.1), the left side of the figure refers to thick-walled rings (the internal radius is indicated), while the right side to thin-walled rings (the mean fiber radius is indicated). Configurations II and IV are the unloaded configurations, corresponding to the stress-free geometry for the thick-walled ring and to the zero transmural pressure geometry in the thin-walled ring, respectively. Notice that in the case of a thick-walled ring the undeformed geometry II differs from the zero transmural pressure geometry III. We refer to Sect. 2.5 for a detailed discussion on the computation of the unloaded configuration. In the last row, we represent generic deformed geometries of the thick and thin-walled ring cross sections, respectively. In the case of the thin-walled ring, we also consider the possibility of section buckling, so that the generic deformed cross section is circular if in pre-buckling conditions (left) or with a non-circular shape if in buckled conditions (right). The same terminology adopted for the radii also applies to the vessel wall thicknesses in the various conditions.
Eq. (9) combined with (8) and (7) and the relative boundary conditions gives with A useful simplification of the expressions in Eq. (10) can be obtained for thin-walled structures, since in this case h u , h 2 u R m,u , with R m,u = (R u + R u,e )/2 (virtual position corresponding to the mean fiber radius). We then obtain the approximations where the second relation represents the well-known Laplace's law. Thin-wall models are considered admissible till γ 1:10 [34], γ being the ratio between the thickness of ring with respect to the radius. As shown in Fig. 5, the venule wall can thus be considered a thin structure, since γ is in the range 1:20 to 1:50. Much different is the situation for arterioles, for which γ 1 : 3, and thus the use of the full expressions in Eq. (10) is required.
In order to derive the expression of the cross section deformed radius, we combine Eqs. (8) with (7) and Eq. (10) (for thick rings) or Eqs. (11) (for thin rings), respectively, obtaining where for thick vessels R denoted the internal radius (blood-vessel interface) while, with a slight abuse of notation, for thin vessels R denotes the mean fiber radius.  [36], data obtained from a meta-analysis of literature studies). Arterioles and venules are designed to withstand different ranges of luminal pressure. The arteriolar wall is thus a thick muscularis layer, while the venular wall is a thin structure. Further, similar, data can be found in [37,38].
The deformed ring thickness can be post-computed from incompressibility, yielding h = R 2 + h 2 u + 2h u R u − R for thick vessels and h = h u R u /R for thin vessels. Functional representation of the Young modulus. Whilst for large blood vessels, especially the carotid, much work has been done, based on experimental measures possibly supported by the use of mathematical models (see, e.g., [3], [39]), there is a substantial paucity of data and models for the representation of the Young modulus in microvessels. Given these premises, if one considers for simplicity E=const, then relations (12) become linear in the pressure, but the corresponding relation transmural pressure vs. cross section area exhibits two non-physiological features: (i) concave form and (ii) absence of saturation at a maximal cross section area for high transmural pressures. It seems then necessary to use a more complex description than a constant, also in view of the different reaction to loads of the components of the vessel wall (collagen, elastin). In this work, we use a linear functional dependence of Young modulus with respect to transmural pressure. This relation is obtained by fitting data from the measurements obtained in [40] by wire myography in small deactivated arteries and veins of the rat mesenteric circulation. The computed steepness of the linear relation is such that the Young modulus passes from a basal value E b at zero transmural pressure to roughly its double when the transmural pressure is increased to 50 mmHg. In the present example of application of the model to the retinal circulation, we set E b = 0.022 MPa for arterioles and E b = 0.066 MPa for venules (basal values chosen as in [9] for the same microcirculatory district). Analogous trends can be obtained also considering different sets of measurements, for example the ones in [39] for human coronary arteries.

Structural model for buckled thin-walled rings
When considering the possibility of reaching buckled configurations, the model must also keep into account the bending actions which actually lead to the loss of axialsymmetry. It is convenient in this context to fix a system of Cartesian axes on the bottom point of the section (see Fig. 6, left). We let s be the arc-length parameter describing the wall mid-line in counterclockwise direction from the origin of the axes and we denote by ϕ = ϕ(s) the angle between the positive direction of the x axis and the tangent to the cross section. The Cartesian coordinates x = x(s), y = y(s) of a point P identified by arc-length s are given by Fig. 6(right) shows an element wedge of arch length ds along with the normal stress σ N , the tangential stress σ T and the bending moment M arising from the pressure loads. According to the hypothesis of thin-walled structure, internal actions have a constant average value in the radial direction. Let K(s) = dϕ ds be the local curvature of the section and K = 1/R u the curvature of the circular undeformed geometry taken as reference configuration. From the approximate theory of curved beams (see, e.g., [41]), the bending moment M = M (s) has the constitutive form where EI is the flexural rigidity, E being the Young modulus (assumed here to be constant and equal to the basal value one for total lack of data) and I = h 3 u /12 the area moment of inertia of the cross section per unit length. The balance of bending moments and forces on the infinitesimal wedge-shaped radial section of ring per unit axial length is given by Combining Eqs. (15) with the Eq. (14) and using the definition of the curvature, yields the nonlinear boundary value system d ds Linear stability analysis of system (16) (see, e.g., [42]) shows that a buckled nonaxisymmetric solution exists for every pressure p t < p t,b , where p t,b = −3EI/R 3 b is the critical transmural pressure corresponding to the lowest energy mode (azimuthal wavenumber equal to 2). When p t = p t,b , the cross section (of radius R b in incipient buckling) loses its circular shape due to physical instability and buckles into an elliptical shape.
For p t < p t,b , progressively, the nearest opposite sides of the section get close, until they touch if the contact pressure p t,c is reached. The contact point becomes a straight line segment in contact if the pressure lowers to the contact line pressure p t,cl . As the pressure is further decreased, the length of the contact line increases and the associated section area tends to zero forming a dumbbell-like shape (see the characteristic shapes reported in Fig. 2).
The buckled configurations have a two-fold symmetry (since they are related to the wavenumber 2), which allows for solving system (16) just in a fourth of the domain. The approach to solve system (16) depends on the value of the transmural pressure, and namely: i) for p t,cl < p t < p t,b , we compute numerically the solution under the hypothesis of isoperimetrical transformations (see also [43] for a similar assumption), by means of the Matlab function bvp4c and using the boundary conditions detailed in [27]. An appropriate choice of the initial guess shape is of fundamental importance to kick in the buckling instability in the computation [19]. We consider a guess shape with a curvature which is a small perturbation (with parameter ε 1) of the curvature of a circle of radius R u , namely K ε = 1 R u 1 + ε cos s R u . This mathematically reproduces the existence in the vessel of imperfections which actually trigger the instability [3]. The value of ε must be tuned accordingly to the imposed transmural pressure in order to obtain a physically coherent solution; ii) for p t < p t,cl , the solution of (16) can be found from that for p t = p t,cl by the similarity transformation [27] ϕ(s) = ϕ cl (s cl ), K(s) = (p t /p t,cl ) 1/3 K cl (s cl ), with the coordinate transformation s = (p t,cl /p t ) 1/3 s cl , where 0 < s cl < s 1 , s 1 being the arc-length of the point of contact in the configuration corresponding to p t = p t,cl .
Once the solution of system (16) has been computed, the non-circular buckled geometry of the section is reconstructed in Cartesian coordinates from Eqs. (13).
The following case study shows an application of the above described model and the hemodynamical consequences of vessel buckling. We consider a single thin-walled vessel (venule) with inlet pressure p in = 40 mmHg and p e = 18 mmHg and we study the flux for outlet pressure p out decreasing monotonically in the range [20,10] mmHg. The vessel has undeformed radius equal to 30 µm and length equal to 370 µm. Simulations are run dividing the vessel into N e = 800 consecutive elements, with progressively smaller elements as the end of the tube is approached. In Fig. 7(left), we show the flux as a function of the outlet pressure. When this latter is decreased, blood flow increases till p out > p e . When p out = p e , the downstream portion of the tube enters into buckling. The flow reaches then a plateau value and it does not depend any more on p out . This trend is in qualitative accordance with the predictions of the Starling resistor model. However, in this latter model only two situations are possible, fully patent or fully closed vessel cross section. The distensible behavior simulated in the present work is more complex, since the conductivities are consistently coupled with the transmural pressure. In Fig. 7(right), we show as an example the 3D configuration of the tube when p out = 10 mmHg. Notice the narrow deformed cross sections in the very downstream portion of the vessel where the low outlet pressure acts. A similar configuration was also observed in [44], where a more complex structural shell model coupled with fluid lubrication theory were used to simulate the experimental setting of the Starling resistor device (notice that in [44] the upstream and downstream cross sections of the tube are maintained fixed).

Computation of the conductivity vs. transmural pressure curve
The knowledge of the geometry of the deformed luminal cross section as a function of the pressure loads allows for computing the element conductivity from relation (3). In detail, we proceed as follows: -if the deformed section remains circular, relations (12) explicitly give the radius of the blood-wall interface. The solution of problem (2) can be then found analytically, and yields the usual parabolic Poiseuille velocity profile [33], from which the conductivity can be straightforwardly computed by integrating (3); -if the section is buckled, problem (2)    Notice how this latter curve significantly differs from the one obtained with the noncircular geometry, especially in the critical area around the onset of the buckling. This motivates us to the explicit computation of the buckled geometry.

Recovery of the unloaded configuration
We conclude the description of the model for a vessel element by dealing with the problem of recovering the unloaded configuration. We assume that the vessel configurations obtained from experimental measurements are circular. Since they do not correspond, in general, to unloaded conditions (configuration II or IV in Fig. 3), we have to solve an inverse problem, whose unknowns are the unloaded configuration itself and the stress field under which the measured deformed configuration is in equilibrium. Let be the generic expression of the structural operator, corresponding to the thick or thinwalled ring models (direct problem). If the undeformed configuration were known, the operator S would compute the measured geometry under given pressure loads. In this context, we have to solve the inverse problem, where the unknowns are the unloaded configuration and the stress field under which the measured deformed configuration is in equilibrium. As, in general, the operator S cannot be analytically inverted, we resort to the fixed-point procedure described in Algorithm 1.
Algorithm 1 : computation of the undeformed geometry Algorithm 1 is similar to the ones proposed in the computational frameworks of [17,45] in biomedical applications, with the introduction in the present case of a relaxation parameter ω r . We have found in our computations that the number of iterations that are actually needed to converge is related to the parameter values, being in particular affected by the wall thickness-to-radius ratio and by the Young modulus basal value.
An example of the application of Algorithm 1 is the following. We start from the deformed geometry (configuration I) of an arteriole with circular cross section of radius R d = 40 µm, thickness h d = 12.8 µm, luminal pressure p d = 40 mmHg and external pressure p e,d = 15 mmHg (data from [28]). The Young modulus is modeled as in Sect. 2.3.1. From Algorithm 1, we obtain the unloaded configuration II represented in Fig. 9(left). Convergence till tolerance 10 −6 is achieved after less than 20 iterations with ω r = 0.3. To give an idea of the importance of reconstructing the unloaded configuration, we also compute configuration III (zero transmural pressure) from II setting p = p e = 15 mmHg and configuration IV, which is the unloaded geometry computed from I using the thin-walled ring model (see again Fig. 9(left)). Observe that configuration IV is only considered for comparison purposes, since the use of the thin-walled ring model is not appropriate with the present value γ = 0.32.
We now apply to configurations I to IV, successively considered as undeformed geometries, the loads p = 20 mmHg and p e = 10 mmHg. In Fig. 9(right), we show the resulting deformation and hoop stress fields.
A significant discrepancy in the stress fields is evident. Configuration I yields stresses which differ of about 1 mmHg with respect to the ones from configuration II (measured vs. unloaded geometry). More significant differences arise if configuration IV is used  Table 2: Geometrical data and stress values for configurations I to IV and percentage variations with respect to values of configuration I. Data refer to the example presented in Sect. 2.5. Row 1 and 3: mean radius and wall thickness of the initial configurations as in Fig. 9(left); rows 2 and 4: percentage variation of the mean radius and wall thickness in the deformed configuration as in Fig. 9(right); row 5: percentage variation of the hoop stress σT at the mean radius of the deformed configuration as in Fig. 9(right). The percentage variation is defined as ∆ r G : instead of II (thin vs. thick structure model). The discrepancies are an increasing function of the magnitude of the transmural pressure, of the individual internal and external pressures, and of the Young modulus (data not reported). Tab. 2 summarizes the geometrical features of each configuration and the percentage difference in results with respect to the ones obtained from configuration I. We conclude this section by noting that in the present work we do not consider the existence of pre-stresses (residualstresses). It is well know that, if cut radially, vessels spring open releasing the residual stress and approaching the zero-stress state which is a circular sector [3]. This aspect is rather delicate and deserves further future analysis.

Extension to a network of microvessels
We now consider the study of the fluid field in a complete compliant microcirculatory network, organized into incoming arterioles, an intermediate capillary bed and draining venules. Referring to Fig. 1 Notice that each element in a vessel has its own cross section area and the number of elements for each vessel may vary but it has its own cross section area. Denoting by Ω i,j f the luminal space of the element E i,j , the fluid domain is given by Junctions between vessel elements and different vessels are simply modeled as nodal points. Let N int be the number of junction nodes. For each node n k , k = 1, . . . , N int , we denote by I k the set of elements which converge in that node. Moreover, we denote by I − k the subset of I k for which n k is the first endpoint of the element, i.e., an element belongs to I − k if its local axis coordinate is such that z = 0 in n k . Similarly, we denote by I + k the subset of I k for which n k is the second endpoint of the element, i.e., an element belongs to I + k if its local axis coordinate is such that z = L in n k . At each node n k , k = 1, . . . , N int , we impose continuity of pressure and conservation of flow (analogue of the electric Kirchhoff's law) At the inlet and outlet nodes n in and n out (physiologically, more than one inlet/outlet can be present in the network, for a total of N bdr boundary nodes), we can apply inlet and outlet pressure values (that is, we impose an overall pressure drop, as in the simulations presented in this work), or an inlet flux and an outlet pressure (or viceversa).
4 Solution procedure

Model summary
The following nonlinear boundary value system of PDEs is to be solved in the compliant domain Ω F : given the connectivity of T , the external pressure, the unloaded configuration and the mechanical properties of the vessels, find the piecewise constant function Q satisfying conditions (19) and the continuous-piecewise linear function p such that in each element it holds dQ dz = 0, where the conductivity σ(p) is determined as summarized in Tab. 1.

Numerical approximation
It is convenient to think that the discrete counterpart of (20) corresponds to the adoption of a primal mixed finite element. In this framework, T represents the "triangulation" of the domain, with elements E. Specifically focusing on the case of prescribed inlet and outlet pressures, we introduce the finite dimensional spaces (see [46] for a similar procedure, albeit in a different context): In Fig. 10, we show the shape function v h,i = v h,i (z) ∈ V h;(g 1 ,g 2 ) relative to node i. This function is linear on each element belonging to I + i ∪ I − i and is such that v h,i = δ i,r , r = 1, 2, . . . , N tot , where we have set N tot = (N int + N bdr ). These characteristics reflect the web-like geometry of the domain.
We let R(p) = 1/σ(p) be the non-negative tube resistance per unit length (observe that σ(p) > 0 in the physiological range) and The numerical solution procedure, including an internal fixed point procedure to solve for the nonlinearities due to the conductivity, is carried out as described in Algorithm 2.
Referring again to Fig. 10 (and omitting for brevity the k superscripts of the internal iteration procedure), we explicitly write the relations obtained from system (23) in Algorithm 2 for a bifurcation with joining node n i , with converging vessel elements E l , E k , E m (here, again for brevity, we have used a shortened notation for vessel elements): Notice that here we have implicitly made use of the fact that functions in V h;(g 1 ,g 2 ) are piecewise linear continuous over T , so that they ensure the automatic satisfaction of the pressure coupling condition. The value of the pressure for all the vessels converging in the node n i is thus uniquely identified by p i . Substituting Eq. (24),(25),(26) (generalized Ohm's laws) in Eq. (27), we end up with a reduced relation in the sole nodal pressure unknowns. Carrying out this procedure for all the nodes, we obtain the linear algebraic system Algorithm 2 : fixed point iteration to compute the fluid-dynamical field on the network h ) on each vessel; compute the cross section geometry from tube law in Fig. 2; update σ where P ∈ R Ntot×1 is the vector of nodal pressure dofs, F ∈ R Ntot×1 is the right-hand side and M ∈ R Ntot×Ntot is the stiffness matrix. Boundary conditions are then enforced in system (28) to ensure uniqueness of the solution. In order to achieve convergence in the internal iteration a relaxation procedure is necessary. We have empirically observed that satisfying a convergence criterion on the pressure but also on the fluxes improves the overall solution. To recover fluxes, we use on each element the the corresponding Ohm's law.
5 Numerical simulations 5.1 A practical instance of "measured" geometry of a microcirculatory network: the case of eye retina vessels The methodology described in the previous sections can address the solution of general unstructured networks. However, in this work we present simulations performed in arteriolar and venular structures with Y-shaped bifurcations and with a mirrored structural organization for arterioles and venules on each side of the central capillary bed. This choice, whilst representing an evident idealization of the real anatomy, allows us to carry out a more systematic discussion on the results, filtering out the effects of the local irregularity and complexity of the geometry.
In particular, we consider asymmetrically branching networks which reproduce the microcirculatory district of the eye retina. Branching is defined according to a modified Murray's Law [47]: letting D f be the diameter of the larger (father) vessel in a bifurcation and D d 1 , D d 2 the diameters of the smaller (daughter) vessels, the following relation is assumed where m = 2.85 is the fractal bifurcation exponent. We assume, as described in [28,29], that the daughter vessels have (possibly different) diameters, given according to where c d 1 /f is a given proportionality coefficient and where D d 2 has been obtained enforcing (29). The generation of the network is continued as long as the vessel diameter is greater than 6 µm. Observe that c d 1 /f = 2 −1/m = 0.784 yields a symmetric dichotomic network with a constant number of branchings leading from the root to each leaf. The length L of each vessel is chosen to be a fractal function of the diameter, according to L = 7.4D 1.15 , as in [28]. Terminal arterioles and venules are connected one-to-one through four parallel capillaries with diameter 6 µm and length 500 µm [28]. Each vessel is then split into equal-sized elements with radius and wall thickness equal to that of the vessel itself. The undeformed geometry is recovered according to the Algorithm 1.
Notice that in this latter configuration the radius of elements of the same vessels may be different due to nonuniform pressure load.
To give an idea of the influence of the asymmetry degree of the network, we consider four different networks with progressively increasing symmetry (that is, with increasing index c d 1 /f ), till reaching a symmetric dichotomic network c d 1 /f = 0.784. In Tab. 3, we report the values of different features of these networks (total number of vessels, min and max route distance of the leaves of the tree, total cross section and equivalent resistance of the network in the measured configuration) for the considered values of c d 1 /f . Shown data refer to the arterial side. The trend of the parameters is due to the increasing homogeneity of the network, which affects the relation between radius and vessel length, and to the constraint of not trespassing the minimum diameter. These elements combined together result into an equivalent resistance which is more than doubled passing from c d 1 /f = 0.5 to c d 1 /f = 0.784.  Table 3: Characteristic values of parameters of networks (arterial side only) generated by different degrees of asymmetry in branching (increasing symmetry moving to the right, 0.784=symmetry). The last column indicates the trend of each parameter for increasing symmetry. Data correspond to a minimal diameter of 6 µm, inlet pressure 40 mmHg, outlet pressure 20 mmHg (values chosen as in [28]).
In Fig. 11 we show an example of network obtained setting c d 1 /f = 0.7. Notice that, whilst the above described fractal networks do not possess, per se, a spatial structure, we endow the network of 3D geometrical coordinates by orienting in the space each daughter branch with respect to the father with elevation and azimuthal solid angles chosen in a range which respects anatomical features. This procedure, on the one side, facilitates the visualization of the network and its physical fields. On the other, more importantly, this allows to locally modify vessel properties or external conditions in a certain 3D spatial region to reproduce physiological and pathological alterations of the baseline values.

Physical and numerical parameters
If not otherwise specified, simulations are run considering coupled arteriolar and venular networks geometrically described as in Sect. 5.1 on each side with c d 1 /f = 0.7. The other parameters are chosen as follows. We set the network inlet pressure to 42 mmHg and the outlet pressure to 18 mmHg, respectively. The inlet arteriole has radius equal to 62µm, the outlet venule equal to 72.5 µm. We set γ equal to 0.32 for arterioles and 0.05 for venules. The first half length of each capillary is considered to belong to the arteriolar network, and thus described as a thick-walled ring, with γ = 0.2. The second half length is considered to belong to the venular network, and thus is described as a thin-walled ring, with γ = 0.08.
Blood viscosity is described according to the model of [48], where it is considered to be a function of plasma viscosity (= 1cP ), blood hematocrit (= 45%) and vessel diameter. In this representation, the viscosity decreases as the vessel diameter decreases till a diameter of about 40 µm, then it starts to increase again with high steepness for smaller vessels. We choose to account for non-circular shapes by using the concept of hydraulic diameter defined as four times the ratio between the vessel cross sectional area and the wetted perimeter of the cross-section [49].
We have studied the influence on the results of the number of element into which each vessel is partitioned. This difference is extremely small when considering the physiological range of typical external and internal blood pressures. It is clear that, were the external pressure be significantly higher than the internal one, using a greater number of elements -eventually with an appropriate adaptive refinement -allows to obtain a more accurate description of the buckling phenomenon. In the following, we set N e = 3 in each vessel.

Simulation results
In the following, we present the results obtained in different test cases which highlight significant network behaviors.
Assessment of the in silico generated networks: comparison with experimental measures. To begin with, we assess the fact that the networks we consider can compute physiologically coherent fluid-dynamical fields. We compare the blood flow we obtain from model simulation considering with the experimental measures in humans performed by different authors in the same diameter range. The external pressure is set to p e = 15 mmHg, corresponding to the intraocular pressure of a healthy subject. The results are shown in Fig. 12, which favorably compares volumetric blood fluxes. Blood velocities (not represented) are also coherent. Pressure field in the network as a function of the interstitial pressure. We computationally solve the mathematical model choosing a uniform interstitial pressure. We set successively p e = [15,18,20,22] mmHg. For ease of presentation, we report separately the results for arteries and veins, even if the simulations are run on the coupled model. We have binned the vessel into classes according to their reference diameters. Class boundaries are obtained by partitioning the range between the minimum and maximum diameter into 20 classes using log10 spacing. Class 1 corresponds to larger vessel, 20 to smaller ones. For each class, we plot the bars whose heights represent the relative number (with respect to the total number of arteries and veins, respectively) of vessels displaying an average pressure in the color-coded range. Colored continuous curves link the height of the bars belonging to the same pressure range. Fig. 13 shows the histogram plots for arteries, Fig. 14 for veins. Observe as, for increasing external pressure, peaks for each pressure range shift towards smaller vessels (larger vessels in the venular network), thus indicating a global increase in blood pressure. Lower pressure ranges progressively tend to disappear. This phenomenon is to be ascribed to the presence of buckled vessels in veins. In the two bottom panels of Fig. 14 we also observe a much more sharp compartmentalization of lower pressure ranges in class 1. Again, this phenomenon is to be ascribed to the presence buckled vessels. This is also the cause of the larger spreading in pressure ranges attained over the diameter classes for increasing p e . These differences are more evident in the step p e =20 to 22 mmHg (p e =18 to 20 mmHg for the venous network) than in the previous steps. Elbow-like turns or "holes" in the enveloping curves are instead due to the asymmetry of the network and the specific bin classification of each vessel.
Network response to local increases of interstitial pressure. We artificially increase the external pressure to p e = 20 mmHg (elsewhere being equal to 15 mmHg) in a delimited portion of the network, located at the interior of a sphere with center in the venous post-capillary zone and radius chosen to include a sufficient number of vessels (gray-shaded region in Fig. 15(bottom row)). This setting simulates, for example, the presence of an edema. In Fig. 15, we plot color-coded net variation ∆p of blood pressure in the arterial and venous networks along with a zoom of the affected regions. In Tab. 4 we report the maximal (with sign) percentage variation along the arteriolar and venular network of blood pressure, blood flux, cross section area and vessel resistance. Albeit no vessel properly enters into buckling instability due to the chosen parameters, a strong increase in resistance is observed in veins due to the combination of area variation and consequent high rise of blood viscosity, strongly diameter depended in that region. Observe also how blood pressure and flow variations are instead similar, this being connected to the continuity conditions enforced at the artero-venous inter-faces. This generates strong flux diversion and "flux stealing" from other vessels [31,50], localized in the four or five generations surrounding the sphere. Observe that, whilst the no arterial vessel is actually included in the sphere, there exists a region of the arterial network which is also affected by perturbations.  Table 4: Maximal (with sign) percentage variations of the arteriolar and venular networks in the specified vessel parameters subsequent to a local increase of external pressure in the spherical region depicted in Fig. 15.

Conclusions
The ability of single body's organs to bring about large selective variations in the rate of blood perfusion relies on the sophisticated regulatory mechanisms of the peripheral circulatory districts. Blood flow regulation is obtained by variation of the vessel diameter, under the effect of both passive and active actions. In this work, we have focused our attention on the first set of mechanisms, investigating the role of geometrical and structural (the so-called "physical") factors in flow regulation. The key point in this process is compliance: elastic microvessels subject to mechanical loads undergo deformations of their wall, altering the shape of the domain offered to blood and, ultimately, resistance to flow. This implies, in turn, a redistribution of the flow itself in the network. The modeling strategy we have proposed in this work uses a simplified description of blood flow and vessel wall-to-flow interaction in order to make feasible computations on large networks. At the same time, important features -peculiar characteristic of these networks-are retained. Blood flow and pressure drop in each vessel element have been connected by a generalized Ohm's law including a conductivity parameter, function of the area and shape of the tube cross section. These latter have been determined in a consistent way by a thin or thick-wall structural model loaded by the internal and external pressure loads. A buckling model is considered in the case of venules, which can experience low/possibly negative values of transmural pressure.
Simulations carried out using the mathematical model show that globally increasing the external pressure causes a global increase in luminal pressure. This process is gradual till buckling occurs somewhere, typically at the outflow, in the venous district. At this point the process has a sort of discontinuity with a much more marked pressure increase. Locally increasing the external pressure, on the other side, has an influence which can be estimated to extend till four or five vessels generations away from the perturbed area. One important point emerges from the above results. Vessels of the same size may experience different intravascular pressure values due to their different location in the network. Hence, even though the vessels embedded in a certain tissue can be classified according to size or branching order, the hemodynamical phenomena which are associated with changes in intravascular pressure should be analyzed and interpreted with caution, especially when considering asymmetrically branching systems.
The present model can be used to investigate alterations in microcirculatory beds due to pathological conditions. Hyperglycemia in diabetes mellitus is known to cause important structural, biochemical and functional changes in the peripheral circulation [51]. Alzheimer's disease, atherosclerosis or vascular dementia can change vascular stiffness. Elevated systolic pressure (hypertension) increases the zero pressure lumen area and the buckling pressure, so that an hypertensive vessel is more likely to collapse than a normal one [52], with many possible hemodynamical consequences. To address these points in a proper way, the following issues should be taken into account in the future: i) inclusion of a more realistic and detailed blood rheology. In our work [30], we have represented blood as a mixture of two fluids, plasma and red blood cells, considering the effects of red blood cell partition at the bifurcations and blood viscosity depending on vessel geometry and variable hematocrit. A similar approach could be extended without too many difficulties to the present model; ii) inclusion of a more realistic model of capillary beds and coupling with the surrounding tissue. If, on the one side, capillaries show a pretty regular mesh-like organization almost ubiquitous in the body, on the other, it is not feasible to face their complete description. Dual mesh techniques as discussed in [53] for cerebral microvasculature can be used to study the dialog between the hemodynamical network and the surrounding tissue. Alternatively, homogenized models for porous media can also be used to effectively face this problem as done in [54] for pulmonary alveoli circulation; iii) inclusion of "acute autoregulation" capabilities in arterioles, which actively respond to stimuli to maintain an adequate blood supply and nutrient delivery to tissue. This topic is naturally founded on having a disposal a model for bloodtissue interaction as described in point ii). Possible regulatory laws could be borrowed from the works by Ursino and co-authors (see, e.g., [55] and references therein) in the case of compartmental models, and from the works by Arciero and co-authors (see, e.g., [56,8]) in the case of a simplified representative 1D/segment model of large/small arteries/veins and capillaries. In these works, the vasoactive response of the arterioles is modeled via "regulatory variables" which depend on myogenic and metabolic stimuli according to phenomenological laws. Similar approaches have been recently adopted to simulate vessel recruitment by David and co-authors (see, e.g., [57] and references therein) for realistic network of the cerebral microcirculation, and by Secomb and co-authors (see, e.g., [15] and references therein).