A physiological model for interpretation of arterial spin labeling reactive hyperemia of calf muscles.

To characterize and interpret arterial spin labeling (ASL) reactive hyperemia of calf muscles for a better understanding of the microcirculation in peripheral arterial disease (PAD), we present a physiological model incorporating oxygen transport, tissue metabolism, and vascular regulation mechanisms. The model demonstrated distinct effects between arterial stenoses and microvascular dysfunction on reactive hyperemia, and indicated a higher sensitivity of 2-minute thigh cuffing to microvascular dysfunction than 5-minute cuffing. The recorded perfusion responses in PAD patients (n = 9) were better differentiated from the normal subjects (n = 7) using the model-based analysis rather than characterization using the apparent peak and time-to-peak of the responses. The analysis results suggested different amounts of microvascular disease within the patient group. Overall, this work demonstrates a novel analysis method and facilitates understanding of the physiology involved in ASL reactive hyperemia. ASL reactive hyperemia with model-based analysis may be used as a noninvasive microvascular assessment in the presence of arterial stenoses, allowing us to look beyond the macrovascular disease in PAD. A subgroup who will have a poor prognosis after revascularization in the patients with critical limb ischemia may be associated with more severe microvascular diseases, which may potentially be identified using ASL reactive hyperemia.

Introduction

Peripheral arterial disease (PAD) refers to the circulatory disorders in the leg as a consequence of peripheral atherosclerosis. Standard diagnostic techniques, such as the ankle-brachial index (ABI) and duplex ultrasound, have performed quite well for the detection of flow-limiting arterial stenoses. However, the influence of PAD is not confined to the conduit arteries but also affects the microvascular function of flow regulation and perfusion reserve [1]. Local perfusion in the tissue bed is dynamically controlled by the vascular smooth muscles of the arterioles, which adjust the wall tension in response to local vasoactive substances. While the standard techniques are poorly suited to assessing the microvascular function, a post-occlusive reactive hyperemia test, as illustrated in Fig 1, has been proposed. In this ischemia-reperfusion paradigm, the blood flow stopped by an inflated cuff would cause hypoxia and accumulation of metabolites in the tissue, leading to functional vasodilation and reduced flow resistance, such that the perfusion after the cuff is deflated can rapidly reach levels several folds higher than the resting condition. Arterial spin labeling (ASL) is a magnetic resonance imaging (MRI) technique that quantifies tissue perfusion in a spatially and temporally resolved fashion. With ASL, patients with PAD commonly demonstrate reduced and delayed reactive hyperemia in their calf muscles following a 5-minute cuffing [2-4].

Fig 1

ASL reactive hyperemia induced by 2 minutes of flow interruption.

TTP = time-to-peak.

However, using ASL reactive hyperemia to detect microvascular disease of PAD remains challenging. First, ASL reactive hyperemia is usually characterized by the peak and time-to-peak (TTP), which are determined based on a single data point of the perfusion response. Because the ASL response has low signal by nature, the noise has great influence on the variability of the two measured indices. Second, the exact vasodilatory mechanisms responsible for leg reactive hyperemia have not been fully understood, which limits the interpretation of abnormal responses observed in PAD. It is not clear if a blunted response implies arterial stenoses, microvascular dysfunction, or a combination of these effects. Lastly, the common 5-minute cuffing protocol was chosen in the past to induce large and prolonged reactive hyperemic arterial flow, for which the presence of flow-limiting stenoses would be highlighted. Given that the microvascular relaxation usually occurs tens of seconds after cuffing, increasing the sensitivity to microvascular dysfunction may require a different protocol.

The goal of this work is to provide the physiological basis of calf reactive hyperemia for data analysis and protocol establishment. We develop a subject-specific model to extract physiologically relevant features in individual responses. Our hypothesis is that a model-based analysis may help explain how the arterial stenoses and microvascular dysfunction affect the overall shape of the ASL reactive hyperemia response in the calf. The specific aims are: 1) to describe the tissue metabolic response during different cuffing durations up to 5 min; 2) to simulate reactive hyperemia under normal conditions and conditions of arterial stenoses and microvascular dysfunction; 3) to fit our acquired reactive hyperemia of cuffing duration-varying experiments in healthy subjects; and 4) to demonstrate the model’s potential of physiological interpretation in preliminary patient studies. If the aims are achieved, this study may lead to a more comprehensive assessment of limb ischemia and better classification of patient groups based on more than the macroscopic lesions in their arteries.

The model

To deal with the challenges in characterization and interpretation of ASL data and standardization of the cuffing protocol, it is desirable to establish a model describing the essential physiological factors involved in reactive hyperemia. Flow models essentially describe the integrated effects of vascular segments connected in series, representing the arteries, microvasculature, and veins. Each vascular segment can be described with a combination of flow resistances (R) and vessel compliance (C). For dynamically regulated flow, the elements of resistance and compliance are made controllable to represent the microvascular responses to factors such as luminal pressure, shear stress, and vasoactive substances. Ursino et al. [5, 6] laid out a theoretical framework to link the physiological factors, including oxygen and metabolites, to the lumped circuit elements for cerebral blood flow regulation in various experimental conditions. However, the modeled reactive hyperemia was far off from experimental reality [7], mainly due to the lack of accurate characterization of the rapid vasodilatory response mechanisms. Adapting this framework, Spronck et al. [8] recently proposed a model suitable for fitting individual cerebral flow responses to experimental maneuvers. Meanwhile, Secomb’s group modeled the range of steady-state skeletal muscle perfusion with extensive study of metabolic flow regulatory mechanisms [9]. In addition, de Mul et al. [10] published a simple model which characterizes PAD patients’ laser Doppler reactive hyperemia with just two time constants, for the rise and decay of reactive hyperemia. Although the model did not take microvascular physiology into account, it demonstrated that a simple RC block is sufficient to characterize the effect of arterial stenoses on perfusion. Therefore, although a physiological model of skeletal muscle reactive hyperemia had not been established, there are plenty of reference materials for use as a starting point. We created the skeletal muscle perfusion model by adapting Spronck’s model and replacing the description of cerebral flow control with ischemia-induced vascular relaxation.

A lumped parameter model of three vascular segments, shown in Fig 2, was used to describe the tibial circulation. Only the arteriolar segment contains adjustable elements, which regulate the flow through the entire modeled system. We solve the lumped parameters of the system as functions of time in order to simulate perfusion-time waveforms. The flow rate in each segment can be derived from Ohm’s law. Specifically, ASL perfusion f measured at the tissue level should correspond to Q in the model: where the parameters are defined in Fig 2 caption. The steady-state baseline flow Q of the system is set to be the group average of our healthy subjects. The segmental pressure distribution of the system at Q is chosen to match the values in the references (see parameters in Table 1 where P was extracted from [11], V from [12], and the rest parameter values from [8]). The governing differential equations to describe the system’s response to flow interruption with time-varying lumped parameters will be derived.

Fig 2

The flow circuit of calf circulation.

P, P and P, denote arterial input, venous output, and intramuscular pressure, respectively. Q and Q denote popliteal arterial inflow and venous output flow, respectively. Q and Q are the flow to the arteriolar and venous circulation, respectively. The switches S1 and S2 represent the flow blocking effect of cuffing. R: popliteal resistance, C: popliteal compliance, P: end arterial pressure for perfusion input to the arterioles; R: arteriolar resistance, C: arteriolar compliance, P mid-arteriole pressure; R: venous resistance, C: venous compliance, P mid-vein pressure.

Table 1

Constant parameters of the model.

Symbol	Value	Unit	Description
Pvo	14	mmHg	Venous output pressure
Pim	10	mmHg	Intramuscular fluid pressure
Qr	144	mL/min	Resting flow rate to the lower leg
Vm	3	L	Estimated volume of one lower leg
ηb	3×10−5	mmHg⋅s	Blood viscosity
ra,0	75	um	Inner radius of unstressed arteriolar wall
ha,0	25	um	Wall thickness of unstressed arteriole
σ1	11.19	mmHg	Constant for elastic tension model
σ2	52.51	mmHg	Constant for elastic tension model
ke	4.5	–	Constant for elastic tension
Tm,0	3	mmHg⋅s	Maximal muscular tension of arterioles
ra,m	128	μm	Radius at maximal muscular tension
ra,t	174	μm	Constant for muscular tension model
nm	1.75	–	Constant for muscular tension model

The lumped parameters

The popliteal artery is modeled as a simple resistor-capacitor unit. The normal arterial resistance is calculated assuming Poiseuille flow: where the arterial radius r = 0.25 cm and length l = 40 cm are estimated from population data, resulting in a resistance of 0.013 mmHg⋅min/mL and an arterial pressure drop P − P = 1.8775 mmHg at Q. The arterial capacitance C is chosen such that the arterial time constant matches the reported experimental values τ = RC ≈ 5 s [10].

Venous resistance is chosen such that a venous pressure drop is P − P = 1 mmHg at Q. With an empirical choice of the venous time constant RC = 10 s, C can be calculated. R and C are taken to be constant [8].

Because the arteriolar circulation is controlled by complex microvascular regulatory mechanisms rather than apparent lumped parameters with fixed values, it is necessary to express R and C with time-varying radius r and consider the relationship between the pressure and radius based on the mechanical properties of the arterioles [6, 8, 13]. In this regard, the mid arteriolar pressure P is related to the wall tension T and radius r by Laplace’s law. The wall tension consists of the passive elastic component T and active muscular component T scaled by the level of activation A. Together these give: where the arteriolar wall thickness h is formulated by: The elastic tension is defined as: where constants are given in Table 1. The muscular tension generated by the fully activated arteriolar smooth muscle at a given radius is defined as: An arteriolar pressure drop ΔP = P − P = 70 mmHg and mid-arteriole pressure P = 50 mmHg are assumed for Q [14]. The resistance R and volume V (related to C) expressed as a function of r is shown in Appendix A. The flow equations of the system.

With the lumped parameters defined, the governing equations of the flow system are expressed as:

Microvascular regulation

The contribution of muscular tension to the total tension is adjusted by its’ activation level A, which is in the range [0, 1] and is defined as a sigmoidal function of the total regulatory influence z:

The regulatory influences are explained in the following. At first, the myogenic and shear stress-based influences were considered because the blood pressure and flow are directly perturbed by thigh cuffing. Myogenic response refers to the reaction of vascular smooth muscle when the vessel is expanded by an increase in blood pressure. Shear stress is the drag force of viscous blood to the vessel wall. However, these two wall-derived mechanisms do not involve vasoactive substances. To date, the mechanisms specific for ischemia-induced vasodilation have not been clearly explained. Fortunately, evidence in the literature reveals that two microvascular regulatory mechanisms are very likely to be tied to ischemia-induced vasodilation. First, it is found that ATP released by red blood cells (RBCs) during hypoxia can activate the purinergic receptors on the endothelial cells and triggers endothelium-dependent vasorelaxation [15]. Second, ischemia-induced imbalance of energetic supply and demand would lead to an accumulation of interstitial adenosine, which can activate the receptors on vascular smooth muscles and induces endothelium-independent vasorelaxation [16]. Therefore, we defined the total regulatory influence as: where i = [myo, sh, ATP, ado] stands for the myogenic, shear stress-based, intravascular ATP-mediated, and interstitial adenosine-mediated regulation, respectively. For each regulatory mechanism, the influence is the regulatory state x multiplied by the respective gain g. The parameter x sets the baseline arteriolar tension, radius, and resistance, which essentially determines the dynamic range for flow increases in the system.

The regulatory state is perturbed by the associated stimulus and varies based on a certain time constant (via a first-order approach): where τ stands for the time constant of each regulatory mechanism. The definition of the stimulus function y is explained in Appendix B. Modeling the regulatory influences. In short, y is the deviation of current arteriolar wall tension from the baseline tension, y the deviation of current shear stress from the baseline shear stress, and so on. Since the wall-derived regulatory mechanisms have been well characterized and are known to be relatively weak when compared with metabolic flow regulation [9], we chose to use fixed values similar to a previous study [8]: g = 1, g = 1, τ = 6 s, and τ = 60 s.

We assumed that intravascular ATP and interstitial adenosine are the main vasoactive substances for metabolic regulations because their production depends on the level of oxygenation in the capillary and tissue, respectively. Therefore, a two-compartment model for the transport of oxygen is required. We adapted Lai’s model and used a constant permeability surface area-product (D) for diffusion of oxygen between the capillary and tissue compartment [17]. The differential equations of compartmental oxygen pressure are: where PO2 and PO2 are the oxygen partial pressure of the capillary and tissue compartment, respectively. O2 and O2 are the oxygen concentration of the inflow arterial and capillary blood, respectively. The solubility of oxygen in blood α is 1.46 μM/mmHg, and volume fraction of capillary v and tissue v are 3% and 97%, respectively. The metabolic function VO2 and the parameter γ = dO2/dPO2, varying with the level of oxygenation, are defined in Appendix B. Modeling the regulatory influences.

The concentration of intravascular ATP C can be calculated by considering the inflow arterial ATP concentration C, oxygen saturation-dependent ATP release rate R, and degradation by the endothelial surface: where C = 0.1 μM, the capillary radius r = 3μm, tube hematocrit h = 0.4, discharge hematocrit h = 0.3, and degradation constant k = 2 × 10−4 cm/s [9]. The ATP release function R is given by in Arciero’s [9]: where R0 = 84 μM/min, R1 = 0.891, P50, = 26.8 mmHg, and n = 2.7.

As ischemia progresses in the tissue compartment during cuffing, the skeletal muscle switches from aerobic metabolism to anaerobic metabolism and the formation of interstitial adenosine C is increased. On the other hand, interstitial adenosine is cleared slowly by cellular uptake back to the tissue with the rate described by Michaelis-Menten equation [18]. Therefore, the dynamic of C is formulated as: where V = 100 μM/min and K = 200 μM are chosen for the adenosine clearance rate [19]. By referring to the concentration in various physiological conditions [20, 21], the adenosine formation rates F in μM/min are assumed here to be: where the critical oxygen pressure PO2 for anaerobic metabolism is assumed to 10 mmHg [22]. The baseline interstitial concentration of adenosine is 0.1 μM.

Model behavior

Simulation of ischemia and reperfusion under normal and diseased conditions is presented in this section to facilitate interpretation of the model behavior. This simulation was developed based on a set of standard parameters reflecting the mean response of the healthy group, as described later in the methods and results sections.

During ischemia, the muscle tissue continues to consume oxygen at the baseline metabolic rate until the PO2 is relatively low, as simulated in Fig 3a. The concentration of intravascular ATP increases with deoxygenation in the capillary compartment, whereas the interstitial adenosine only starts to increase when PO2 is reached and anaerobic metabolism begins, as delineated in Fig 3b. The influences of each regulatory mechanism are drawn in Fig 3c. While the total regulatory influence z continues to be enhanced by the accumulation of metabolites, the vascular smooth muscle is almost de-activated after 2 minutes of ischemia, as shown by the level of activation A in Fig 3d.

Fig 3

Ischemic responses of the tissue and arterioles.

Modeling the responses during 5 minutes of ischemia in the muscle, including (a) the oxygen tension, (b) concentrations of the vasoactive substances, (c) the influences of vasoregulation mechanisms, and (d) the summed influence z and level of activation A of vascular smooth muscle.

The modeled reactive hyperemia is shown in Fig 4. Because of the rapid nature of intravascular ATP-mediated vasorelaxation, an evident response can be induced with just one minute of ischemia. As illustrated in Fig 4a, increasing ischemic duration would slightly enhance the magnitude of the modeled response, but ischemia longer than 2 minutes would mostly prolong the modeled reactive hyperemia. Arterial stenoses simulated by doubling R would reduce the response magnitude regardless of ischemic duration, as depicted in Fig 4b, 4c and 4d, because it reduces the input perfusion pressure to the microcirculation. Doubled R would also double the arterial time constant, resulting in longer rise time. Microvascular dysfunction does not delay the rise time, on the other hand. Depending on the ischemic duration used, microvascular dysfunction does not necessarily reduce the magnitude of the responses to a large extent. As simulated by reducing g in the model, microvascular dysfunction appears to cause earlier attenuation, suggesting that the dysfunctional arterioles tend to constrict. When comparing microvascular dysfunction alone with normal function, or combined micro- and macro-vascular disease with stenoses alone, the differences are more evident in the responses to 2-min ischemia than longer ischemia. This result of modeling may be explained by the fact that as the concentration of interstitial adenosine increases in longer ischemia the vascular smooth muscles eventually relax, regardless of the endothelium-dependent ATP signal. Therefore, the model suggests that a protocol using shorter ischemic duration may be more sensitive to endothelium-related microvascular dysfunction.

Fig 4

The effects of ischemic duration, arterial stenosis, and microvascular dysfunction on reactive hyperemia.

(a) Normal responses to various ischemic durations and diseased responses to (b) 2-min, (c) 3-min, and (d) 5-min of ischemia. Arterial stenosis was simulated by doubling the popliteal resistance R, and microvascular dysfunction was simulated by reducing g to 40%. “Combined” includes both effects at the same time.

Materials and methods

Subjects and data acquisition

The perfusion data was acquired previously in a separate work of sequence optimization, and a portion of it was summarized in S1 Supplementary Material. In brief, the study included the mid-calf ASL reactive hyperemia of 7 healthy subjects (age <30 yrs) and 9 patients with PAD. All of the healthy subjects underwent 2, 3, and 5 minutes of arterial occlusion via cuff inflation (220–250 mmHg); five of them also participated in a 1-min arterial occlusion experiment. Each of the patients only underwent 2-min cuffing for the more affect leg, which was determined based on self-reported symptoms and physicians’ dictation. The median age of the patient group was 74.5 yrs (51 to 83). The ABI ranged from 0.3 to 0.87. Seven patients had claudication, and the remaining two had rest pain.

The ASL technique was flow-sensitive alternating inversion recovery (FAIR) with single-slice steady-state free precession for acquisition, performed at 3T (MR 750 by GE, Milwaukee, Wisconsin) with an 8-channel cardiac receive array placed in the calf region. Sequence parameters were: field of view = 16 cm × 16 cm, matrix size = 64×64, slice thickness = 1 cm, single-slice steady-state free precession (SSFP) acquisition, flip angle = 70°, post-labeling delay (PLD) = 1.4 s, bolus duration TI1 = 700 ms. Image acquisition was repeated every 3.35 s. The images were reconstructed and processed in MATLAB. Regions of interest were manually drawn to enclose the entire muscle area excluding the large vessels. The difference signal ΔM(t) was converted to perfusion f(t) as follows: where M0 is the signal of fully relaxed blood, T1 the blood longitudinal relaxation time, and the blood volume partition coefficient λ = 0.9mL/g. Reactive hyperemia was characterized by the peak and TTP.

Modeling and fitting

The flow regulation system described previously, including the lumped parameters and oxygen transport model, was implemented using MATLAB (see S1 Model File, the entire model called RHmodel). Modeled reactive hyperemia was generated by solving the governing differential eqs (7)–(9) and (12), using the built-in solver function ode15s. Subject-specific model parameters were identified as described in the next two paragraphs. The values of these parameters were estimated by searching for the minimum least-squares fit between the modeled and ASL-measured perfusion, using the built-in optimization function lsqnonlin. The quality of fit was assessed by calculating the reduced chi-square: where σ2 is the variance of baseline perfusion, N the data length, and n the number of free parameters.

For the healthy group, [f, g, τ, g, τ, τ] were chosen as the free parameters. The logic lies in the approach to address the individual differences in the resting perfusion f (= Q/V) and reactive hyperemia. First, we attributed the variation of resting perfusion in healthy individuals to the baseline arteriolar tone, reflecting the variation in capillary permeability and tissue metabolism, not the variation in resistances of the arteries or veins. Therefore, f was varied to determine x for the baseline arteriolar tone. The group average resting perfusion was aligned with a positive x = 0.45 chosen empirically. Next, we assumed that the differences in reactive hyperemia result from the variation in the sensitivities and time constants of ATP and adenosine in addition to the arterial time constant (reflecting the difference in arterial compliance if the resistance is fixed). Therefore, the steps of fitting the cuffing duration-varying reactive hyperemia are: 1) use [f, g, τ, τ] with fixed g and τ to fit the responses to 3- and 5-min cuffing; and 2) use [f, g, τ] with the resulting g, τ, and τ from step 1 to fit the responses to 1- and 2-min cuffing. The parameters g and τ were not included as free parameters in step 1 because their influences on the response to long cuffing were extremely small and a wide range of their values can be used, making the resulting values unstable and not meaningful. Likewise, only the more sensitive parameters for short cuffing were included in step 2. The allowed range of free parameters in the fitting process are: 2 ≤ f ≤ 10 (mL/100g/min), 1 ≤ g ≤ 20, 1 ≤ g ≤ 15, 6 ≤ τ ≤ 24 (s), 12 ≤ τ ≤ 60 (s), 1.5 ≤ τ ≤ 30 (s).

For the patient group, [f, g, τ, R, τ] were chosen as the free parameters. R was included as a free parameter because these PAD patients have various degrees of arterial stenoses affecting reactive hyperemia. Adenosine parameters were not included as free parameters because of their minimal influence on perfusion response to 2-min cuffing. Their values were fixed to the average values of the healthy group. The allowed range of free parameters in the fitting process are: 3 ≤ f ≤ 8 (mL/100g/min), 1 ≤ g ≤ 12, 12 ≤ τ ≤ 36 (s), 1.5 ≤ τ ≤ 25 (s), and 1 ≤ R ≤ 6 (starting from here R is expressed as the ratio to the calculated standard value from Eq (2)). The ranges of these parameters were narrower than those used in the healthy group because the physiological ranges had been known after fitting the healthy responses.

Statistical analysis

The peak, TTP, and model parameters that resulted in the best fit to the responses were compared between the healthy subjects and patients using an unpaired t-test. A p-value of 0.05 was assumed to indicate statistical significance.

Results

The model could generate ischemic-duration-dependent responses very similar to the measured reactive hyperemia in healthy subjects. The model achieved reasonable fits with one set of parameters for all responses to various durations of cuffing in each subject, as shown in Fig 5. The quality of fit, measured by , was affected by baseline drift, sudden jumps, and small oscillations of the ASL signal in some cases. Noise level varied between recordings. Looking through the healthy group, the model did not necessarily fit better to the large and long responses (5-min cuffing) than the transient responses.

Fig 5

Fitting the responses of a healthy subject with the model.

The resulting parameters are: [f, g, g, τ, τ, τ] = [5.59, 14.67, 3.49, 18.0, 44.5, 3.05].

The parameter values of the healthy group were summarized as means ± standard deviations: f = 4.95 ± 0.92 mL/100g/min, g = 15.66 ± 2.05, g = 4.18 ± 2.41, τ = 15.86 ± 2.66 s, τ = 43.20 ± 21.11 s, and τ = 4.25 ± 1.16 s. The parameter values and of each fit were reported in S1 Dataset. The parameters related to adenosine had a broader range than other parameters, with coefficients of variation (CV) 27.4% and 49.6% for g and τ, respectively. The rest of parameters were relatively consistent within the group (CV<20%).

The model fitting of diseased responses to 2-min cuffing was shown in Fig 6. Eight out of the 9 responses had ranged between 0.59 to 1.39, indicating reasonable fits; response 3 had a of 4.17, indicating a poor fit (see S1 Model File for each fitting). Each recording had a different level of noise. Response 3 was very small and only three data points were two standard deviations above the baseline perfusion, which was not sufficient to determine the five free model parameters.

Fig 6

Fitting of patients’ responses to 2-min cuffing with the model-generated response.

The fitting with exceptionally large error is plotted as the blue line.

The model parameters were more variable within the patient group than the healthy group. As compared in Fig 7, responses in the patient group were significantly lower in peak perfusion, but not in TTP (p = 0.102), than those in the healthy group. With the model-based analysis through fitting, the patient group showed higher R and lower g than the healthy group. The p-values were 0.285, 0.098, and 0.461 for comparison of f, τ, and τ.

Fig 7

Characteristics of reactive hyperemia induced by 2-min cuffing.

Left: characterization using apparent indices. Right: characterization using model-derived parameters. * denotes statistical significance (P < 0.05).

Discussion

Establishment of the model

We sought to provide an explanation to the dependence of reactive hyperemic characteristics on ischemic duration using the model. Although reactive hyperemia is often related to endothelial function, which is commonly assessed by flow-mediated dilatation via the effect of shear stress, it is unlikely that shear stress plays a significant role in ischemia-induced vasodilation. The first reason is that flow interruption should lead to shear-dependent vasoconstriction, not vasodilation. Second, the change in shear stress induced by cuffing is independent of the ischemic duration and hence should not enhance the perfusion response if ischemic duration increases. The same arguments can be applied to the pressure-related myogenic regulatory mechanism. Indeed, it has been suggested that these wall-derived mechanisms should be relatively weak and work against the metabolic flow regulation [9]. To identify the mechanisms more relevant to ischemia, we searched for mechanisms closely tied to the dynamics of the deoxygenation process during ischemia. Studies of magnetic resonance spectroscopy have shown that myoglobin desaturation begins at 100–120 s and plateaus at 300–420 s in the course of acute ischemia [23, 24]. The initial delay of myoglobin desaturation in the tissue pool is due to the contribution of oxygen from the blood pool. At the later stage of ischemia, high level of myoglobin desaturation corresponded to the onset of anaerobic metabolism indicated by the change of phosphocreatine in the tissue. These two-compartment deoxygenation dynamics were incorporated in our model (Fig 3a) to control the metabolite-derived regulatory mechanisms.

In the capillary compartment, we linked vascular tension to the oxygen saturation-dependent release of ATP by RBCs. The theory of RBCs being the oxygen sensor for local flow regulation has long been established [15]. Our simulated concentration of ATP was in the submicromolar range, consistent with intravascular measurements in human hypoxia experiments [25]. The estimated time constant of ATP-mediated regulation τ was 15.4 ± 2.7 s, which is very close to the reported values (8–16 s) from direct observation of exposed arterioles in an animal study [26]. The function of the ATP receptor on the endothelium was found to be attenuated in type 2 diabetes [27]. This finding corresponded to a reduction in g in our model and was simulated as microvascular dysfunction. Consistent with the simulation, reduced g was found in the patient group.

In the tissue compartment, the concentration of interstitial adenosine was modeled. Increased formation of adenosine during hypoxia results from the elevated activity of ecto 5’-nucleotidase converting muscle-released AMP into adenosine [28]. However, it is difficult to quantify the exact reaction rate and concentrations of AMP and adenosine because a wide range of values was reported, potentially because real-time interstitial measurements in vivo remain challenging [29]. Hence, we only roughly estimated the adenosine formation rate and described it as a simple 2-mode function in Eq (18). The high variability of g and τ in the fitting of normal responses (Fig 6) may result from such simplification of adenosine formation. Furthermore, other metabolites, such as lactic acid and hydrogen ions, may also be involved in the tissue compartment for extended ischemia. Therefore, the current model may be less accurate in explaining/simulating the responses to longer ischemia.

With intravascular ATP and interstitial adenosine incorporated, the model clearly demonstrated ischemic duration-dependent responses. The rapid magnitude increase for responses to short ischemia is very likely to be contributed by the ATP released by deoxygenated RBCs. On the other hand, longer ischemia appeared to prolong the responses by accumulating adenosine in the tissue. The relationship of the peak and TTP on ischemic duration in the modeled reactive hyperemia agreed with the previously measured characteristics in the healthy group. From a statistical perspective, however, some of the acquired responses did not entirely agree with the model. ASL data may contain not just metabolic perfusion response but also imaging noise, motion artifacts introduced by cuff inflation and deflation, and higher order vessel properties causing flow oscillation during recovery. These details were not considered in the model, which may be the reason for poor fitting in some cases. Previously, Spronck’s model was used to fit individual cerebral flow response averaged by repeated recordings, which presumably decreased the experimental variability associated with these influences. Therefore, our model may fit better to averaged individual responses, which should be more representative of metabolic flow responses. The fitting may be very likely to fail when the diseased response is very small. However, cases with such small responses may essentially indicate severe microvascular disease and model analysis is not required.

Major findings

An important contribution of this work is that the model delineates the different effects of arterial stenoses versus microvascular dysfunction on reactive hyperemia. Arterial stenoses can limit response peak to a great extent. However, there are already many tools to assess arterial lesions. The use of ASL reactive hyperemia was motivated by the hope to look at microvascular function. According to our simulation, a shorter ischemic duration can better detect microvascular dysfunction. This finding reflects the often overlooked fact that the microvasculature typically reacts quickly to ischemia and a short ischemic protocol can be used to measure the quickness of reactivity. To date, 5-min cuffing is recommended for reactive hyperemia measurement due to its lower variability in the peak and TTP when compared with shorter cuffing. From a signal perspective, our simulation agrees with this viewpoint that the 2-min cuffing response may be too short to provide enough data points to measure the exact peak difference caused by arterial lesions. However, our model simulation suggested that the peak and TTP are more useful for detecting the influence of macroscopic lesions rather than quantifying microvascular dysfunction. The model clearly indicated that the unique potential of ASL reactive hyperemia for microvascular assessment is highlighted when using 2-min cuffing. Therefore, we proposed to use a 2-min ischemic protocol for more sensitive microvascular assessment.

With 2-min cuffing, the PAD group demonstrated lower peaks and slightly longer TTP than the healthy group (Fig 7), which is similar to the previous findings using 5-min cuffing [4]. A significant difference in TTP may be found if subject size increases. Another strategy to raise significance is to reduce variability by taking model fitted peak and TTP, since the raw peak and TTP is based on a single data point so is highly susceptible to noise. Fitting the responses to a known mathematical function, such as a gamma-variate function, has been demonstrated for reactive hyperemic blood oxygen-level dependent imaging [30]. While the function may also fit ASL reactive hyperemia of the healthy subjects, the diseased responses may not resemble a gamma-variate function, as exemplified by response 9 in Fig 6. More importantly, an empirical function does not contain physiological meaning, let alone explaining the dependence of reactive hyperemia on cuffing duration. With our model-derived analysis, we found that some of the model parameters not only clearly differentiate the responses of the two groups but may also suggest direct pathophysiological implications. First, the increased R may represent arterial stenoses. Second, the arterial time constant τ is a function of both resistance and compliance. Since the two factors are not necessarily tightly correlated and have the opposite effect on τ as they progress with PAD, a wide range of values in patients is not surprising. Lastly, the g, representing the microvascular endothelial sensitivity to ATP in our model, was significantly lower overall in patients but also had greater variation within the patient group. This result suggests that different levels of microvascular disease exist within the patient group and may be assessed by ASL reactive hyperemia with our model-based analysis. The variation in the microvascular parameter may have critical clinical implications.

It is known that a large portion of the patients with critical limb ischemia does not achieve limb salvage 1 year after revascularization, which could be caused by more advanced microvascular diseases. This work strongly demonstrates the potential of using ASL reactive hyperemia for microvascular assessment in such a patient population, which may be useful to identify the subgroup who requires more effective treatment for their microvascular disease.

The scope of this work was determined to focus on the mathematics and physiological concepts about reactive hyperemia. As a first step, we only explored the perfusion characteristics in healthy subjects and patients with known disease. Testing the clinical relevance of the derived perfusion indices will be the next focus. The perfusion measure for the arterial resistance can be tested by a comparison between the more and less affected leg of the same patient. A bilateral comparison can exclude the variability caused by the differences in subjects’ cardiac output, blood pressure, blood viscosity, and other systematic factors, and the results may be more meaningful than a comparison between patients and age-matched non-PAD subjects. On the other hand, a standard microvascular assessment that can be used to validate the perfusion measures of microvascular function does not yet exist. However, the relevance of the microvascular measures to the outcome can be tested by a post-revascularization prognostic study.

So far, the model only explained the data from small subject groups. With more data accumulated, the patients can be stratified based on their clinical indications, and the perfusion indices between subgroups can be compared. Also, the dependence among model-derived parameters and the correlation between the parameters and clinical indices may be further examined. The model description of reactive hyperemia should not be generalized for exercise-induced active hyperemia, as capillary recruitment in exercise can greatly change the permeability surface area-product. Since the heterogeneity between calf muscle groups was not the focus here, perfusion of the whole mid-calf cross section was used for the establishment of the model. Using perfusion data from smaller ROIs with higher noise in ASL does not help extract consistent physiology for a general description of the response to ischemia. Although the model fitted reasonably well to our data, the quality of ASL data remains critical. Large SNR and minimum motion artifact can reduce bias in the fitting. The experiments of 2-min thigh cuffing were well tolerated by patients, but repeated recording could be affected by previous ischemia, especially in patients with reduced washout rate of metabolites.

Conclusion

We established a detailed physiological model to provide a description of the flow regulatory mechanisms involved in ASL reactive hyperemia of calf skeletal muscles. The model demonstrated that combining the dynamics of intravascular ATP released by RBCs and interstitial adenosine can explain the evolution of response characteristics with increasing durations of cuffing in healthy subjects. Model simulation suggested influences of microvascular dysfunction distinct from those of arterial stenoses on reactive hyperemia. Microvascular dysfunction may be detected more easily in a 2-min cuffing protocol rather than longer cuffing, whereas the effect of arterial stenoses may be more consistently observed using 5-min cuffing. Model-based analysis of the patient data showed that R and g were significantly different between the patient and healthy group. While increased R in patients is normally expected, the variation of g within the patient group may suggest that different amounts of microvascular disease existed. Therefore, we propose next to test the utility of ASL reactive hyperemia and our model-based analysis in the patients with critical limb ischemia who require assessment of their microcirculation in addition to their arterial lesions.

Selected symbols

f: Resting perfusion; g: Gain of intravascular ATP-mediated regulation; g: Gain of interstitial adenosine-mediated regulation; R: Popliteal artery resistance; τ: Response time to adenosine; τ: Response time to ATP; τ: Popliteal artery inflow time constant; : Reduced chi-square.

Appendix

A. The flow equations of the system

The arteriolar resistance R and volume V are defined as: where the constants K and K are chosen such that at r = r (baseline), R = ΔP/Q, RC = 1 s, and V = CΔP.

The arterial flow eq (7) is derived from the definition of compliance being a change of volume due to pressure change: where V is the volume of the popliteal artery. The venous flow eq (9) is derived likewise.

The arteriolar flow eq (8) is derived by taking the time derivative of Eq (23) and expressing the volume change as the difference between the input Q and output flow Q:

B. Modeling the regulatory influences

The stimulus functions of wall-derived regulatory mechanisms are taken from Spronck’s model [8]. The myogenic stimulus is calculated based on the deviation of arteriolar wall tension T from the baseline T0, multiplied by a scaling factors s: where s = 3.33 (mmHg⋅cm)−1.

The shear stress is proportional to Q/r3. Therefore, the shear stress stimulus is defined as: where s is chosen such that at rest the baseline flow and arteriolar radius will result in zero stimulus.

Likewise, the stimulus functions of metabolite-derived regulatory mechanisms, including the intravascular ATP and extravascular adenosine, are defined as the deviation from the baseline concentrations C and C multiplied by scaling factors: where s and s are empirically chosen to be 2.5 and 0.61 μM−1, respectively.

The concentrations of oxygen in the capillary O2 and tissue compartment O2 are the sums of freely dissolved oxygen and oxygen bound to hemoglobin and myoglobin, respectively: where C = 2.33mM is the concentration of hemoglobin, C = 365μM the concentration of myoglobin, w = 80.95% the fraction of water in blood, and w = 78% the fraction of water in the tissue. The parameters γ = dO2/dPo2 for Eqs (13) and (14) are:

With O2 defined, the skeletal muscle metabolic function is calculated: where K = 0.7 μM, and V is calculated from the individual baseline perfusion.

The change of intravascular ATP concentration Eq (15) is derived by considering the total amount of ATP in a single capillary with inflow, release and degradation: where the capillary length l, flow rate of a single capillary q. Since capillary flow can be derived from perfusion and capillary volume fraction, this equation can be simplified into Eq (15).

RHmodel and perfusion data to produce Figs 3, 4, 5 and 6.

This zip-file contains functions and perfusion data in MATLAB’s format.

(ZIP)

Click here for additional data file.

Result values of model fitting in the healthy group.

(XLSX)

Click here for additional data file.

Result values of direct characterization and model fitting in Fig 7.

(XLSX)

Click here for additional data file.

Acquisition of ASL reactive hyperemia and result images.

This supplementary file includes ASL sequence settings, image processing, quantification of perfusion, result perfusion maps, and comparison of perfusion-time curves between the healthy and patient group.

(PDF)

Click here for additional data file.