Introduction

The microcirculation is the essential ‘end point’ of the cardiovascular system, consisting of a vast network of microvessels perfusing the body’s organ tissues, whose main function is delivering oxygen and nutrients and removing waste products. At the local level, the microcirculation continuously regulates the levels of blood flow and pressure across its network to satisfy metabolic demands, to redistribute hydraulic loads and to promote inflammatory processes. To this end, a sophisticated hierarchical control system (intrinsic, metabolic and neurohormonal) operates by regulating the diameter of the microvessels, causing vaso-dilation or vaso-constriction, which in turn modulates their hydrodynamic resistance. When a significant increase in tissue metabolic demand occurs, a large fraction of the microcirculation needs to be recruited for coordinating the vaso-dilation. To achieve this, the vaso-dilation arising from the lower microcirculation can ‘ascend’ into feed arteries by means of an electrical signal conducted along the endothelial cells (ECs), hyper-polarizing vascular smooth muscle cells (SMCs) causing relaxation.11

Since microcirculation represents the ‘mesoscale’ functionally bridging the systemic circulation to perfused tissues, microvascular dysfunction (MVD) ultimately compromises oxygen delivery to the end organ(s).121 Indeed, it is within the microcirculation that the earliest signs of several cardiovascular diseases manifest themselves. For example, the Firefighters and Their Endothelium (FATE) study8 showed that microvascular health represents a powerful independent predictor of cardiovascular events in primary prevention. Thus, assessment of microvascular function may not only have utility for providing novel insights into the pathophysiology of the patient, but also presents an important opportunity for early disease detection and risk stratification.7

Microvascular function can be evaluated invasively at the level of the end-organ using vaso-active agents combined with pressure/flow wires and non-invasively using advanced imaging techniques such as computed tomography (CT), magnetic resonance imaging (MRI) and positron emission tomography (PET).22,108 Invasive approaches can indeed selectively partition macro- from microvascular function but are not suitable nor feasible in all individuals. While non-invasive CT, PET and MRI address many of these concerns, exposing participants to ionizing radiation, confined spaces, and/or strong magnetic fields also limits widespread adoption. As such, there has been considerable focus on the development of safe and cost-effective alternative methods for measuring microvascular reactivity.123 Accordingly, much attention has been placed on the peripheral vasculature, as it is easily accessible and reasonably reflects vascular function in other end-organs of interest.6,21,90 Although not yet established in clinical practice, methods such as reactive hyperemia provide prognostically-significant indicators, capable of differentiating clinical phenotypes.106,107 In its simplest form, reactive hyperemia represents the magnitude of limb reperfusion following a brief period of ischemia induced by arterial occlusion (see Fig. 1).

Figure 1
figure 1

Microvascular conditions before, during and after occlusion. The size of the arrow indicates the magnitude of blood flow. The vasodilation in resistance arteries may have a smaller entity than in arterioles due to the limited distance that can be reached by the ascending hyperpolarizing signal.

Full size image

Fundamentally, the approach measures microvascular vasoreactivity to metabolic substances produced in response to tissue ischemia. Here we consider reactive hyperemia to be limited only to the initial increase in flow that immediately follows the restoration of flow. Indeed, the blood flow response after this initial period is affected by regulatory mechanisms that fall outside of pure ‘microvascular function’. Multiple methods exist for evaluating reactive hyperemia, including limb distension by venous occlusion plethysmography, blood velocity/flow via Doppler ultrasound of an upstream conduit vessel, kinetic changes in tissue oxygenation by near-infrared spectroscopy (NIRS),12,38,50,66,105 and tissue perfusion by near-infrared diffuse correlation spectroscopy (NIRS-DCS)13; each of which with its own strengths and weaknesses.107 As detailed in this review article, the determinants of reactive hyperemia are multi-factorial and difficult to quantify. Thus, the best measurement approach is one that accounts for the greatest number of such complexities.

Computer modeling can be used as a complementary tool for hypothesis testing, making predictions and quantifying the dynamics underlying vascular regulation that cannot be directly observed during in vivo experimentation. Here, we report on various modeling techniques used to describe reactive hyperemia dynamics across the circulation, and describe the effects of various assumptions on hierarchical blood flow control system. The review concludes with an integration of knowledge, relating the theoretical modeling with existing data collected by our group and others (by way of peripheral ischemia-reperfusion), identifying key opportunities for future discovery.

Haemodynamics

Occlusion of a conduit blood vessel, like the brachial artery, has a direct negative impact on the resulting pressure and flow distribution downstream of the occlusion. The reduction in pressure quickly propagates from the occlusion site down to the level of the exchange vessels, ultimately impairing oxygen delivery. The ensuing tissue ischemia results in microvascular vasodilation (intended to correct the error signal), such that when the occlusion is reversed, the ensuing blood flow is markedly elevated (i.e. reactive hyperemia).

Flow in (Large) Compliant Vessels

The exchange of forces between blood and vascular wall and its resulting displacement can be evaluated by employing detailed three-dimensional (3D) fluid-structure interaction models53 or simplified one-dimensional (1D) blood flow modeling approaches which consider field variations only along the main flow (axial) direction.88 The latter methodology is less accurate (especially around localized anatomical details), but has many competitive computational advantages and can be easily adopted for vessel networks.5 This makes 1D blood flow modeling the best option for describing the fluid mechanics in reactive hyperemia, which involves propagating phenomena along the vasculature with time scale of many seconds. Furthermore, through 1D blood flow modeling it is also possible to perform wave intensity analysis which allows the quantification of pressure waveforms travelling forward to the microcirculation and backward to the heart.88 In arteries (and arterioles) the blood behaviour can be generally approximated as homogeneous and Newtonian since the size of the red blood cells carried by the plasma is considerably smaller (10 times) than the vessel diameter. Furthermore, flow is also generally considered incompressible and laminar, with a Poiseuille velocity profile. In 1D blood flow modeling, each compliant vessel can be treated as axisymmetric, with blood velocity (u), pressure (P) and flow described as continuous variables along its axial direction z. These quantities are averaged over the cross-sectional area (A) and their variation along the radial direction is considered negligible. The Navier-Stokes equations for 1D blood flow in compliant vessels can be expressed in terms of cross-sectional area and velocity averaged over the cross-section:

$$\begin{aligned} \begin{aligned} \frac{\partial A}{\partial t}+\frac{\partial (A u)}{\partial z} = 0, \\ \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial z}+\frac{1}{ \rho } \frac{\partial P}{\partial z} +\frac{\mu u}{\rho }\frac{8\pi }{ A} = 0, \end{aligned} \end{aligned}$$
(1)

where t is time, \(\mu\) is the fluid dynamic viscosity, \(\rho\) is the fluid density, while \(Q=Au\) is the (volumetric) flow rate. It is worth noting that (1) can also be written in terms of flow rate and pressure26 or cross-sectional area and flow rate.112 The mechanics underlying the vascular wall deformation appears to be complex, mainly due to vessel visco-elastic properties and the capacity to produce active tone for diameter regulation. To describe the interaction between blood and vessel wall, different approaches can be used.31,89,113 In the simplest case, the fluid pressure is related to the cross-section of the vessel via a linear function with respect to the luminal diameter (D=\(\sqrt{4A/\pi }\))

$$\begin{aligned} P=P_{\text{ext}}+\beta (\sqrt{A}-\sqrt{A_0}), \end{aligned}$$
(2)

where \(P_{\text{ext}}\) is the external pressure from the surrounding tissue, \(\beta\) is a parameter representing the wall elasticity and \(A_0\) is the unstressed cross-section area. However, modeling reactive hyperemia requires a vessel wall model able to describe the hyperpolarization-induced dilation and then the resulting (compliant) structural response to hyperemic flow. Given such complexity, wall mechanics models derived from conservation laws retaining mechanobiological features33,86 are preferable over tube laws which are purely phenomenological; indeed the latter, for this specific application, would require the introduction of several non-physical parameters. It is also noted that, if desired, wall viscoelasticity can be integrated into the blood pressure-wall deformation law by using more complex constitutive models3,15,96 with a consequent decrease in computational efficiency. The haemodynamic features (viscosity and density) and the geometric (diameter and length) and structural (stiffness) wall properties can be made vessel specific and can reflect different physiological and pathological conditions such as ageing or hypertension.34,88 The blood flow variables described by system (1) and (2) can be computed by employing a broad variety of numerical schemes including finite differences, finite elements and finite volumes and can be extended to large vessel networks by imposing mass and momentum conservation at the interface between vessels.26,82,88,118

Biphasic Flow in Microvessels

During reactive hyperemia, flow in the microvessels is altered from its physiological range due first to the sudden pressure reduction induced by the upstream occlusion and then by arteriolar luminal expansion consequent to the wall relaxation driven by ischemic tissues. Due to its morphology and function, the downstream vasculature constitutes the site of major blood pressure drop along the cardiovascular system. The work by Secomb111 provides a detailed characterization of the flow through microcirculatory networks. In these microvessels the Reynolds numbers are < 1, and therefore the blood flow can be described with a good level of accuracy as incompressible Stokes fluid, for which the convective component is neglected. Blood flow can still be described by (1) but most assume a rigid microvessel wall, which implies considering only the momentum conservation equation for relating blood flow and pressure. Furthermore, the biphasic nature of flow requires a modeling framework that accounts for the main rheological properties of its components. Since they are concentrated in the central part of the tube, RBCs travel faster than the surrounding plasma (referred as Fahraeus effect) and the shorter RBC’s transit time with respect to the total blood flow has a profound effect on oxygen transport and exchange with tissues. The discharge hematocrit (\(H_{\text{D}}\)), which is the ratio between RBCs volume flux and total blood volume flux along the tube, appears higher than the tube hematocrit (\(H_{\text{T}}\)), which is the volume occupied by RBCs over the tube volume at a certain time instant (see Fig. 2).

Figure 2
figure 2

Hematocrit in blood vessel. The discharge hematocrit, being the volume fraction of RBCs in blood delivered through the vessel by the flow, is obtained as ratio the vessel volumetric rate of RBCs and volumetric blood flow (\(H_{\text{D}}\)=\(\dot{V}_{\text{RBC}}/\dot{V}_{\text{blood}}\)). The tube hematocrit is given as ratio between the volume occupied by RBCs and total blood volume in the vessel (\(H_{\text{T}}\)=\(n_{\text{RBC}}{V}_{\text{RBC}}/{V}_{\text{blood}}\), where \(n_{\text{RBC}}\) is the number of RBCs in the vessel).

Full size image

The distribution of RBCs along a microvascular network can be computed by means of the advection equation.25 The mean velocity of RBCs (\(u_{\text{RBC}}\)) can be estimated if the ratio between \(H_{\text{T}}\) and \(H_{\text{D}}\) is known (since \({H_{\text{T}}}/{H_{\text{D}}}\)=\(u/{u_{\text{RBC}}}\)), for which Pries et al.99 identified the following correlation:

$$\begin{aligned} \frac{H_{\text{T}}}{{H}_{\text{D}}}=H_{\text{D}}+(1-H_{\text{D}})(1+1.7 \text {exp}(-0.35D)-0.6 \text {exp}(-0.01D)). \end{aligned}$$
(3)

The latter expression was then modified by Pries et al.102 to account for the reduction in diameter due to the endothelial surface layer. In tubes with diameters below 300 \(\mu\)m, the viscosity of blood reduces with decreasing diameter. This phenomenon, discovered by Fahraeus and Lindquist49 (referred as Fahraeus-Lindquist effect) appears to be associated with the tendency of RBCs to migrate away from the vessel wall.111 The resulting blood viscosity can be evaluated as a product between the plasma viscosity and the relative apparent viscosity (\(\mu _{\text{rel}}\)), with the latter accounting for the effect of the hematocrit. Pries and co-workers defined a set of correlations which relate apparent viscosity with discharge hematocrit, and diameter of microvessel for in vitro100 and in vivo conditions.99,101 In the latter case, the relative apparent viscosity is defined as:

$$\begin{aligned} \mu _{{{\text{rel}}}} = \left[ {1 + (\mu _{{45}}-1 )\frac{{(1 - H_{\text{D}} )^{{C_{\mu } }} - 1}}{{(1 - 0.45)^{{C_{\mu } }} - 1}}\left( {\frac{D}{{D - 1.1}}} \right)^{2} } \right]\left( {\frac{D}{{D - 1.1}}} \right)^{2} , \end{aligned}$$
(4)

in which \(\mu _{45}\) and \(C_{\mu }\) are experimentally derived correlations dependent on the vessel diameter. The presence of cell-free layers within the microvessels is also responsible for the phase separation effect, which consists in a non-proportional partition of RBCs and plasma flows at diverging bifurcations. The experimentally-derived correlations by Pries et al.98,102 represent the reference approach for including this phenomenon within vascular networks. An alternative RBCs phase separation model was subsequently introduced by Gould and Linninger,61 although its benefits over the one by Pries et al.98,102 have been debated.103

Opposite to the continuum approach, blood flow across microvessels can also be described via motion of multiple discrete particles51,56 and the effect associated to the particulate nature can be assessed for different channel sizes.10 To transit across the smallest capillaries, the RBC membrane undergoes significant deformation, with the intensity of applied forces on the capillary lumen playing a significant mechanobiological role.64 We expect that the surge in blood flow during reactive hyperemia substantially changes the RBC-capillary lumen interaction. Variation in RBC number and velocity across the capillary bed can be experimentally analyzed by using intravital video recordings76 but the resulting pressure field remains uncertain. Discrete multi-scale RBC modeling approaches1,28,52 can be used for reproducing the mechanical and rheological properties of the cell and their interaction with plasma. Casquero et al.27 used a NURBS-based fluid-structure interaction framework for assessing the deformation of RBCs, which were modeled as compound capsules (or solid membrane) where the nucleus is treated as bulky deformable solid. These multidimensional models allow the direct assessment of the physics underlying the interaction between RBCs and microvessel wall but they require efficient parallel algorithms62,95 for being scalable to realistic and physiologically relevant problems.

Lumped Circulation Models

Lumped or zero-dimensional (0D) modeling represents an approach in which the system’s quantities are not (directly) space-dependent but depend only on time. Their governing equations can be directly derived from formulations of high-fidelity/space-dependent models (such as full 3D heart chambers, arteries and veins) and the space-dependent characteristics are encapsulated through a set of parameters which are generally identified using experimental data.119 0D models may be employed alone for describing the full cardiovascular system or parts of it and can be represented through electrical analogy. Lumped models can also be used in combination with higher-fidelity counterparts to represent those vascular territories where a low level of space details is required.

Compartmental lumped models20,78,92 split the vasculature into different flow compartments, each of them characterized by specific parameters representing vascular features such as flow resistance and compliance. In this case, the full haemodynamic profiles across the vasculature can be obtained by simply computing a set of (time-dependent) ordinary differential equations. This approach involves very low computational costs and allows an easy integration and/or comparison with quantities averaged over space obtained from some experimental imaging techniques (such as doppler ultrasound or NIRS-DCS). Therefore, due to its simplicity, 0D modeling should be considered as first-step for modeling blood flow variables during reactive hyperemia. On the other hand, lumped modeling alone cannot adequately capture the wave-propagation phenomena underlying the reactive hyperemia response and therefore does not provide a suitable platform for a full mechanistic understanding of the different space-dependent interactions between vascular territories. Nevertheless we re-iterate that, for various questions related to reactive hyperemia, lumped modeling should be considered for the description of the whole system’s dynamics, with space-dependent models introduced as ‘sub-models’ for further refinement if necessary.

Vascular Network Definition

The topology of the fluid network can either be reconstructed from images17,76,115 or generated by using branching (either deterministic or stochastic) algorithms.2,74 From CT or MRI scans it is possible to obtain, after segmentation, networks of large blood vessels17 but diameters smaller than 500 \(\mu\)m are generally below the standard angiography spatial resolution. On the other hand, microvascular networks can be digitally reconstructed, down to the capillary scale, by leveraging recent developments in high-resolution microscopy.115 To facilitate the integration of experimental datasets in network models, methods were developed for classifying microvessels without requiring any flow information115 or estimating flow rate with unknown boundary conditions55 The microvascular arrangement and consequent levels of blood perfusion strongly depend on local morphological and functional properties of the tissue such as \(\hbox {O}_2\) consumption rate, level of venous drainage and presence of anastomosis.131 Skeletal muscles are perfused by bundles of parallel capillaries which transport blood from a common terminal arteriole to the same post-capillary venule. Mendelson et al.76 used intravital video microscopy to show that capillary modules may present substantial structural differences in terms of topology, as well as functional heterogeneities in RBCs haemodynamics within the same module and across modules. This study proposed a new paradigm for capillary flow organization which is pivotal for explaining the complex flow regulation occurring during skeletal muscle contraction. Capillary modules are indeed interconnected to form continuous columns, defined as capillary fascicle, that align with the structure of the muscle fascicle.76 Blood flow through the arterial-venous vasculature (arteries \(\rightarrow\) arterioles \(\rightarrow\) capillaries \(\rightarrow\) venules\(\rightarrow\) veins) is complex, dictated by varying vessel wall morphology and geometry. Due to this, multi-scale network models can be used to accurately capture the blood flow transition across the mesoscale (large to microvessel).2,74

To compute haemodynamic variables along a portion of the vasculature, it is necessary to define the boundary conditions at the network extremities to represent the remaining part of the cardiovascular system that is not explicitly modelled in space (see for instance4). For computing blood flow in large arteries during reactive hyperemia, the inlet boundary condition could be imposed either at the heart level or at a (proximal) brachial artery site (as long as the inlet is not beyond the occlusion site). However, the latter case is practically less feasible than the former as it requires the prescription of both flow and pressure (which are difficult to measure simultaneously) for not omitting the potential contribution of backwards reflections. The left ventricle (or the proximal aortic node) can therefore be used as the inflow prescription point,16,54,88,113 with its complex pumping function simplified using a lumped model, which is able to capture the essential features of inflow data.

To close the fluid network, it is always necessary to define appropriate outflow conditions, which represent the haemodynamics characterizing the downstream systems. Venous hydrodynamic conditions (i.e., position of the limb) may play a significant role in the experimental settings,13 as they directly affect the downstream pressure level. Therefore, inclusion of a part of the venous system into the computer model may be useful for reflecting changes in the functional state of its components. Others have also proposed using a closed loop flow model, in which the arterial network system is coupled with a venous counterpart at the heart and capillary levels.83,89 While this holistic approach incorporates additional features, like the venous network, translation of this model to a patient-specific context is hindered by the large numbers of parameters which need to be identified/measured. Due to their simplicity and adaptability, lumped Windkessel models represent the most popular approach for modeling vascular territories beyond the main system of interest.89 These lumped models can be combined in series to form complex interconnected compartments, allowing to establish a link between their functional states and the outflow conditions. Alternatively, more complex methods can be used to represent the downstream circulation, including the structured (or fractal) tree approach,91 its derived versions,37,94 and the porous-media based approach.32,36 While these models encapsulate more details of the downstream network, they come with higher computational costs.

Coordinated Regulation

Tissue ischemia generates hyperpolarization of the smallest microvessels (capillaries mainly) which is then conducted upstream through EC gap junctions to relax the arteriolar vasculature and favour a more consistent/systemic increase in blood flow. In the time period between tissue ischemia generation and the full recovery of the system, intrinsic pressure and flow control mechanisms simultaneously work together to restore homeostasis; the entirety of which is extremely complex and difficult to quantify.24,73

Generation of the Hyperpolarizing Signal

Tissue metabolism involves the production and consumption of different substances, some of which have been traditionally associated with hypoxia (i.e., adenosine, \(\hbox {K}^+\), \(\hbox {CO}_2\), lactate and \(\hbox {H}^+\)). In reactive hyperemia, the mismatch between oxygen demand and supply is caused by a sudden decrease in oxygenated blood flow and not by a change in tissue metabolic activity, which instead drives functional hyperemia. The review by Rosenberry and Nelson107 reports an exhaustive description of the cellular pathways involved from the stimulus generation to vessel dilation. Accordingly, the stimulus originating the hyperpolarization is yet to be defined but most likely involves a combination of \(\hbox {K}^+\), ATP, bradykin, \(\hbox {H}_2\)\(\hbox {O}_2\) and epoxyeicosatrienoic acids. These molecules cause vessel hyperpolarization by activating different vascular channels including inwardly rectifying potassium channels and \(\hbox {Na}^+\)/\(\hbox {K}^+\) pumps.38 Recent findings75 on the capillary EC function showed that membrane hyperpolarization may trigger a complex cascade of \(\hbox {Ca}^{2+}\) signaling pathways which include cross-talk between intracellular \(\hbox {Ca}^{2+}\) stores and transient receptor potential channels. In a hypoxic environment, RBCs may also contribute to microvascular dilation by releasing substances into the bloodstream such as ATP.46,47,58,72,117,133

In the computer model, the originating stimulus can be represented either as a single variable in a set of kinetic equations describing oxygen/tissue signaling or simply as perturbation of the concentration of species in the extracellular space. In the reference model by Moshkforoush et al.81 the effect of cerebral metabolic activity on capillary ECs was represented as an increase in extracellular K\(^+\) concentration, which in turn affects transmembrane K\(^+\) channels and membrane potential. Hyperpolarization of capillaries can be modeled by employing electrochemical models35,81,114,126 describing the dynamics of intracellular species (such as \(\hbox {Ca}^{2+}\), K\(^+\) and Na\(^+\)) and membrane potential in stimulated ECs. These models are considered as ‘space-homogeneous‘ (or lumped/0D) because their variables represent quantities averaged over the cell volume and therefore only dependent on time. In the seminal work by Wiesner et al.126 the effect of an exogenous agent (thrombin) was simulated by combining the EC \(\hbox {Ca}^{2+}\) dynamics with kinetics models representing the time evolution of ligands involved in the receptor-activation pathway. Modeling studies that quantitatively characterize the function of a component of this orchestra, such as,35 are therefore instrumental for deciphering the emergent cell dynamics features. The capabilities of space-homogeneous models can be further extended by adding transport equations for describing diffusion of species127 or models representing cytosolic micro-domains80 for assessing the interaction between sub-cellular components. Space-homogenous models allow an easy extension to cluster of cells for computing the resulting intra-cellular signaling. On the other hand, this would be hard to achieve in case a multi-dimensional approach for cell dynamics was used. Therefore, cellular space-homogeneous models represent the most trivial option as bottom level of a multi-scale vessel framework. Furthermore, \(\hbox {O}_2\)-induced release of ATP from RBCs across capillary networks can also be accounted for in the model by following the work by Goldman et al.57,60 However, a recent theoretical study104 suggested that, with respect to other O\(_2\)-driven regulatory pathways, the role of vasoactive molecules release by RBCs is limited.

Conducted Vasodilation

Ascending vasodilation allows to export the local capillary hyperpolarization to a widespread fraction of the microvasculature. The endothelium serves as the predominant cellular pathway for signal conduction as the hyperpolarization is transmitted from cell to cell along the vessel wall through gap junctions.11 According to Welsh et al., the conduction of the signal is enabled by a defined pattern of charge movement along the vascular wall which is the result of the interplay between tissue structure, gap junction resistivity and ion channel activity.125 Despite the vessel axis being the main direction of propagation, the hyperpolarizing signal is also transmitted into surrounding SMCs through myoendothelial junctions to promote vascular relaxation.122 Endothelial \(\hbox {Ca}^{2+}\)-activated K\(^+\) (SKCa/IKCa) channels seem to play a key role in tuning electrical conduction along microvessels by modulating signal dissipation through changes in transmembrane resistance.14 Experiments in rat cremaster arterioles indicate that vasoactive conduction is favoured by circulation of autocrine and paracrine mediators ATP and K\(^+\), which is increased under hyperemic conditions.44

Early efforts in modeling electrical communication across vascular branching include the works by Diep et al.,42 in which vascular intercellular gap junctions in the skeletal muscle were represented by ohmic resistors. For investigating the signal transduction involved in conducted vaso-reactivity, Kapela et al. integrated their previous modeling efforts69,114 into a multicellular model of a rat mesenteric arteriole which comprises a detailed description of ECs-SMCs units, coupled via non-selective gap junctions.70 This framework was used in a subsequent study71 for estimating from experimental measurements the electrical resistance associated to cell membrane and gap junctions. More recently, the same group presented an experimentally-validated computational framework for reproducing the propagation of the hyperpolarizing signal across a cerebral microvascular network.81 This provides some key information regarding how the activity of different ion channels can impact on the ascending signal and its re-generation. However, this framework neither quantifies the dynamic changes in microvascular compliance nor indicates whether the metabolic demand is satisfied but it represents a solid modeling base to start from. Alternatively, Arciero et al.9 developed a multi-compartmental framework based on a previous autoregulation wall mechanics model to represent the adjustments in diameter across a microvascular network induced by vaso-conducted metabolic response. In this work the stimulation considered was the hypoxic-induced ATP release by RBCs and the hemodynamic model was coupled with Krogh-type models for representing the mass exchange with tissue.

Vessel Dilation

The conducted changes in membrane potential spread from capillaries to arterioles. Indeed, in rabbit skeletal muscle18,19 and studies on brain function63 capillaries are able to control local blood through pericytes contraction. In the study by Horn et al.,65 a feed artery occlusion caused a heterogeneous flow distribution in capillary beds with the presence of capillary no-reflow, and this appeared as result of the complex interaction between myogenic and metabolic mechanisms occurring at the arteriolar level. In simulated capillary modules, RBCs flow appeared to be independent from the network resistance.77 This indicates that RBCs flow across capillaries is mainly regulated at control sites located at-pre and post-capillary levels. Overall, the role of single capillaries on skeletal muscle local flow regulation still needs to be fully elucidated but its extent seems minor with respect to the upstream microcirculation.97 Experimental evidence indicated that also venules and veins are endowed with a contractile apparatus, and therefore they may play a role in the pressure re-distribution occurring during and after the arterial occlusion.39,43,120 However, their high blood volume capacity suggests a minor role of these vessels especially in the first part of haemodynamic response after occlusion release.

The level of dilation of each microvessel depends on the magnitude of the hyperpolarization signal that has reached its location and by the ion channels expression levels of its wall constituents.87 A change in the wall membrane potential promotes the increase in vessel diameter in potentially two ways. First, as direct effect, hyperpolarization of the SMC causes a decrease in intracellular \(\hbox {Ca}^{2+}\) which ultimately diminishes the activity of myosin light chain kinase with reduction of cross-bridges formation. A second effect on the contractile machinery may be mediated through ECs’ \(\hbox {Ca}^{2+}\) dynamics. Indeed, in these cells, the decrease in membrane potential causes an increase in intracellular \(\hbox {Ca}^{2+}\) which promotes the formation of autocoids such as NO which ultimately affect the cross-bridges by phosphorylating the myosin light chain phosphatase. Both pathways, by reducing the fraction of attached actin and myosin filaments, leads to distension of the SMC and consequent vessel dilation. However, Crecelius et al.38 showed that, in the forearm, there is no combined role of NO and PGs in peak RH. Therefore, EC \(\hbox {Ca}^{2+}\) dynamics appears to be crucial for the conduction of the hyperpolarizing signal along the endothelium and to the adjacent SMCs while its paracrine function is less relevant. In skeletal muscle resistance arteries, \(\hbox {Ca}^{2+}\) sensitization pathways are also known to play a role in the myogenic response,79 but these seem to be present only above a certain blood pressure level (\(\sim\) 60 mmHg).93

To accurately model the change in vessel diameter due to membrane hyperpolarization a mechanistic multi-scale model is required. The contractile machinery model needs to translate the biochemical information from SMC \(\hbox {Ca}^{2+}\) signaling into tissue deformation. For this purpose space-homogenous models for describing cellular signaling seems the most viable choice. The SMC \(\hbox {Ca}^{2+}\) dynamics model developed by Kapela et al.,68 based on 26 partial differential equations, provides a comprehensive description of the interaction between several ionic channels upon different agonist stimulation. Pioneering work by Yang et al.129,130 bridges the multi-scale components of SMC contractility. This framework, devised for computing the myogenic response in rat cerebrovascular arteries, includes a description for the cross-bridges kinetics and corresponding cell length variation, which establishes a mechanistic link between intracellular \(\hbox {Ca}^{2+}\) variation to structural vessel deformation. Similarly, Layton’s group developed other reference models29,45 able to accurately mimick the myogenic response in renal afferent arterioles. On the contrary, Carlson and Secomb23 presented a minimalistic vessel wall mechanics model that, without including any cellular/sub-cellular component, provides a good agreement for both passive and active tension against myograph experimental data from single isolated microvessels. This type of approach is extremely attractive due to its simplicity and limited computational cost but a dynamic interplay with ECs network carrying the hyperpolarizing signal remains challenging to establish. The multi-scale models following the approach by Murtada et al.33,84,85,86 provide a rigorous way to encapsulate the information regarding the actin-myosin interaction into a tissue/continuum mechanics level by using an active strain energy function. This class of models represents the best option for coupling the resulting wall deformation and stress with the blood flow within the vessel.

Regulated Networks and Clinical Data Integration

To date, there is a limited number of models able to describe the coordinated haemodynamic response to arterial occlusion through a multi-scale vascular system. The work by Yamazaki and Kamiyama128 is one of the earliest attempts to simulate flow-mediated dilation in a conduit vessel by including multiple scales (from cell to vessel) into one modeling framework. Jin et al.,67 extended this work by including a realistic blood flow network model. This model closely approximates experimental recordings made at the conduit artery level but its representation of the downstream microvascular dynamics remains limited. Some studies59,109,110 have also assessed the impact of microvascular diameter changes on blood rheology and tissue perfusion across different vasculatures. It is however important to note that, despite such models allowing quantification of flow characteristics upon imposed variations in microvessel resistance, there is no direct link with a realistic coordinated regulatory response.

Due to their intrinsic simplicity, compartmental lumped models represent the standard approach for evaluating the blood flow dynamics across different levels of the microvasculature, without directly describing propagation phenomena. For each compartment, vascular changes due to regulation can be represented by varying (in time) its lumped parameters (such as resistor and capacitor, see Fig. 3) according to more or less complex internal variable evolution equations.30,40,116,124,132

Figure 3
figure 3

Example of regulated blood flow compartmental model with potential clinical data integration. \(P_{\text{ai}}\): inlet arterial pressure; \(Q_{\text{a}}\): arterial flow; \(R_{\text{a}}\): arterial resistance; \(C_{\text{a}}\): arterial compliance; \(Q_{\text{ma}}\): micro-arterial/arteriolar flow; \(R_{\text{ma}}\): micro-arterial/arteriolar resistance; \(C_{\text{ma}}\): micro-arterial/arteriolar compliance; \(Q_{\text{c}}\): capillary flow; \(R_{\text{c}}\): capillary resistance; \(P_{\text{vi}}\): inlet venous pressure; \(P_{\text{ext}}\): external pressure from surrounding tissue. \(R_{\text{ma}}\) and \(C_{\text{ma}}\) are modulated to reflect vascular changes due to regulatory mechanisms (green arrows).

Full size image

de Mul et al.40,41 were among the first to describe the fluid-dynamics occurring during post-arterial occlusion reactive hyperemia by using compartmental modeling, through which vascular territories are split into arterial and capillary systems. In Reference 40 measurements from Laser–Doppler perfusion monitoring were used for estimating two time constants characterizing the dynamics, which differed across different types of patients (healthy vs unhealthy). Likewise, Vo et al.124 modelled reactive hyperemia response to venous occlusion and by using NIRS data. Solovyev et al.116 defined an hybrid model, combining a lumped model for blood flow with an agent based model of skin injury, to assess the role of blood flow on the development of pressure ulcers in spinal-cord injury patients. Recently, a multi-compartmental model incorporating oxygen transport, tissue metabolism, and vascular regulation mechanisms was proposed by Chen and Wright30 to characterize and interpret MRI-derived microvascular reactive hyperemia (by arterial spin labeling) in calf muscles of human volunteers. Through all these modeling approaches it is generally possible to evaluate quantities that can be directly associated with in vivo experimental measurements carried out in human patients. However, this type of approach requires the condensation of most regulatory mechanisms within few model parameters, which might not have a clear physical meaning and that in most of the times need to be estimated by fitting against experimental data. At the same time, some of these compacted model parameters can be used for a straightforward discrimination of different types of cardiovascular profiles. As such, given the complexity of reactive hyperemia, this class of models is likely most appropriate for capturing the interaction between vascular compartments and their contribution to the regulated haemodynamic response.

Discussion and Concluding Remarks

There is a pressing need to enhance our understanding on microvascular function and its regulation. Advances in this direction will ultimately lead to the development of new diagnostic, progress monitoring and therapeutical approaches for diseases characterized by (micro)vascular dysfunction. Through this work, we sought to provide a guide for modeling reactive hyperemia, which represents a standard method for assessing microvascular reactivity. Based on the current knowledge on reactive hyperemia, we reported a series of models that can be used for describing physiological phenomena that occurs at different space and time scales, spanning from cellular ionic currents to systemic circulation flow dynamics. Several issues hinder the accurate modeling of reactive hyperemia. First and foremost, from an experimental point of view, isolating single vascular control mechanisms and vascular compartment contribution to the resulting haemodynamics is technically very challenging if not impossible. On the other hand, it is possible to measure the combined response at locations across different experimental conditions, in order to best inform the computer model. The use of multiple simultaneous recordings of blood flow via Doppler ultrasound as well as RBCs flux via NIRS-DCS may also be helpful in the identification of model parameters and quantification of their uncertainties.48

There is need to define, test and use suitable modeling approaches that, despite the complexity of the system analysed, are capable of capturing the relevant interactions between vascular compartments. Furthermore, there is need for computational models that can efficiently deal with different scales in space and time as this aspect is essential for the successful translation of the computer model into clinical practice. Connections between components of the hierarchical multi-scale models play a crucial role. Further investigation and validation is needed before modeling may inform clinical decision making. To this end, high-fidelity models that leverage large existing datasets augmented by machine learning/data-driven approaches may be needed.