Modeling of the hemodynamics in the feet of patients with peripheral artery disease

: To simulate the hemodynamic effects in the feet in response to a thigh cuff occlusion, we have developed a multi-compartmental model in which the circulatory system for the leg is represented by its electrical equivalents. Dynamic vascular optical tomographic imaging data previously obtained from 20 patients with peripheral artery disease (PAD) and 20 healthy subjects is used to test the model. Analyzing the clinical data with the support of the model yields diagnostic specificity and sensitivity in the 90-95% range, significantly higher than previously reported.


Introduction
Peripheral artery disease (PAD) affects approximately 12 million people in the US and is considered the leading cause of low-extremity amputations, responsible for over 70,000 amputations annually [1].PAD is caused by an accumulation of plaque in the vessels of the lower extremities, which leads to stenosis and reduction in blood flow.If not properly treated, PAD can lead to critical limb ischemia with consequent cells death, ulcers, and gangrene [2,3].The disease affects mostly the elderly, diabetics, and patients with kidney problems.Early symptoms include claudication (pain while walking), leg numbness, or weakness and coldness in the lower leg or foot [4].
Different techniques are currently used to diagnose PAD.The most commonly employed method consists of determining the ratio of the blood pressure measurements taken in the lower leg and the arm.This so-called ankle-brachial index (ABI) can be obtained noninvasively at a low-cost.However, it has shown to have lower sensitivity in patients with arterial calcifications, such as diabetics (about 20-40% of all PAD patients) and subjects with renal disease [5,6].Another non-invasive technique widely used is duplex ultrasound, which can detect vessel narrowing in the upper thigh, but does not provide information about the oxygenation saturation of tissue in the foot, and small vessels are difficult to visualize [7,8].A third method is angiography, which can provide a map of the vascular tree in the legs.But this technique either relies on nephrotoxic contrast agents (when magnetic resonance imaging is employed) or radiation exposure (when X-ray computed tomography is used), which limits its use in patients with renal failure or for a continuous monitoring of the condition.Moreover, none of these methods provides information about the distal perfusion in the foot, where most of the ulcers occur [7][8][9].
To overcome some of these limitations, we have developed a vascular optical tomography imaging (VOTI) system, which measures transmitted and reflected light intensities and allows to estimate hemoglobin distribution inside the feet [10,11].In a previous clinical study [12], we demonstrated that dynamic time-dependent image features can be used to distinguish healthy from PAD patients with accuracy, sensitivity, and specificity of over 80%, depending on the diagnostic cut-off points.
To further improve our ability to correctly diagnose PAD and to better understand the factors that affect dynamic optical measurements during a dynamic thigh cuff occlusion, we present here a mathematical model of the behavior of the hemoglobin in the lower leg during a thigh cuff occlusion.Employing this model to previously collected clinical data from 20 PAD patients and 20 healthy subjects, provides new insights into the physiological underpinnings of the observed effects.Furthermore, we show that diagnostic sensitivity and specificity can be improved.

Dynamic model of blood flow in foot
To assess the dynamic accumulation of blood in the foot of healthy and PAD patients, we have adopted a multi-compartmental derivation of the Windkessel models.In general, these models (also called zero-dimensional or lumped models) simulate the blood behavior in the vasculature with an electrical circuit analogy [13][14][15][16][17][18][19][20][21][22].Their main advantages are that they rely on simple first-order differential equations [17,20], reduce the number of parameters to consider [19], and can be easily related to physiological parameters.Therefore, the equations encountered are relatively easy to solve and help to understand the global distribution of pressure, flow, and blood volume under various physical conditions [18,20,21].
Employing this approach, the vascular tree can be represented by adding in series or parallel basic lumped-parameters compartments, each characterized by its resistance (which represent the resistance to the blood flow), capacitance (which represent the vessels' elasticity and compliance) and inductance (which represent the blood inertia) [17,18,21].Consequently, the blood pressure and flow are modeled with a voltage and a current respectively.The relationship between the pressure in each vascular compartment and the blood flow through it is represented using Ohm's Law (V = R x I) and a set of first-order differential equations represents the capacitors (C x dV/dt = I) and inductors (L x dI/dt = V).
The arteries, arterioles, capillaries, venules and veins are then characterized by different components and values.The arteries are elastic or viscoelastic and the blood flow is pulsatile.Therefore, arteries display all typical effects related to resistance, capacitance, and inductance [18].In the arterioles and capillaries, the resistance dominates over the other parameters [18].In the venules and veins, the compliance is much higher than in the arteries [13,22] and the blood flow is relatively steady; therefore, an RC combination is usually considered sufficient to describe hemodynamic in the vein compartment [18].
However, since we are interested in the accumulation of blood in the foot over a minutelong cuff inflation, fast changes (on a sub-second or second-scale) in the blood flow only play a minor role.We can, therefore, simplify the model by removing the inductors, which are used to represent fast changes in the blood flow [18].This reduces the total number of components present in the system and we are left with four RC sub-sections: 1-large and small arteries; 2-arterioles; 3-capillary bed and venules; and, 4-veins.Each resistor is connected in series to the following resistor, whereas the capacitor is connected in parallel to the following compartments, forming a ladder circuit as depicted in Fig. 1.Previously published models focused on estimating the shape of the blood pulse or comparing two different system states [14][15][16][17][18][19].These studies did not consider the effects of a dynamic venous occlusion in the leg caused by the application of an external force.To model the effects of this venous occlusion, we changed the existing models by introducing two timevariable resistances for the main veins and arteries sections to simulate the inflation and deflation of the thigh cuff.Considering that the cuff inflation causes a reduction of the blood flow through the blood vessel at the cuff location, one will observe an increase in the pressure before that location.Consequently, the resistance will increase during an inflation.After the deflation, this pressure decreases again to the resting-phase value and hence one will observe a decrease in resistance.
Furthermore, the cuff was operated manually, thus maximum pressures of 60 mmHg and 120 mmHg were not obtained immediately.To model this effect, we explored linear and exponential functions to simulate the increasing/decreasing values for the resistances.We found that exponential fits provided a better match with the data.
Finally, we considered the possibility of a partial delayed arterial occlusion.It is well known that, in a cuff experiment, veins are compressed before arteries become compressed.We, therefore, employ a delay parameter accounting for this fact in our simulations.Taking the above considerations into account, the following system of four first-order differential equations was obtained: Here R denotes the resistances: R A for the small and large arteries, R AO for the arterioles, R C for the capillaries, R VN for the venules, and R V for the veins; C denotes the capacitances: C A for the large arteries, C AO for the arterioles, C C for the capillaries and venules, and C V for the veins.V denotes the voltage drops on the capacitances as indicated by the subscripts; R V (t) and R A (t) denote the dependency of the veins and arteries to the tight cuff occlusion.

( ) ( ) ( )
(5) In these equations t start_i is the starting time for the cuff inflation, t end_i is the ending time of the cuff inflation, t start_d is the starting time for the cuff deflation, t end_d is the ending time of the cuff deflation, Δt is a delay before the arteries are affected by the cuff and R CUFF-Vm represent the additional vein resistance value when the cuff is completely inflated.
R CUFF-V and R CUFF-A are the resistances in the veins and arteries generated by inflating and subsequently deflating the cuff.We assume here that R CUFF-V and R CUFF-A are correlated and X is a multiplication factor between 0 and 1 that links the cuff effect on the veins and arteries: if X is equal to 0 it means that the cuff inflation is not affecting the arterial system, if X is equal to 1 it means that its effect on veins and arteries is the same.As veins are located more towards the skin and arteries are found deeper inside the tissue, a 60mmHg cuff will mostly affect the veins (X is closer to 0).Arterial occlusion typically occurs at pressures above 100mmHg (X closer to 1).An example of the time-dependent functions R V (t) (red curve) and R A (t) (blue curve) is shown in Fig. 2. In this example X = 0.2, meaning that the cuff induces a 5-times larger resistance in the veins than in the arteries.This system of Eqs.(1) through 5) can be solved by employing the Euler method [21] with time-step Δt smaller than the one used to characterize the heartbeat.The model was implemented in MATLAB.A graphic user interface was created to facilitate the manipulation of the different parameters.The heartbeat (V Heart (t) in Eq. ( 1) was modelled using typical electrocardiogram and physiological pressures [13,21,23,24]: heart rate of 75 beats/min, systolic pressure of 120 mmHg and diastolic pressure of 80 mmHg.The cuff inflation and deflation times used were set to mean values, respectively 15 seconds (t end_i -t start_i ) and 5 seconds (t end_d -t start_d ), obtained from experimental tests.The arterial delay time Δt was set to 2/3 of the time necessary to inflate the cuff, since low pressure values do not affect the arteries and the pressure during the cuff inflation follows the exponential behavior shown in Eq. (2-5) and Fig. 2.
The value range to be used for the main resistances and capacitances of the lumped model is not uniquely defined in literature, and many different values have been used and measured in different situations, while no one has been proven superior [18,22].After a literature review, comprising all the references cited in this paper, Table 1 shows the range of values we considered for resistances and capacitances.

Experimental data used for testing the model
To test our model, we used clinical data previously obtained from 20 healthy volunteers and 20 PAD patients.The design and results of that study were previously reported in reference [12].Here we provide a brief summary.
The VOTI system used during the experiments employed a total of 34 optical fibers, 14 of which were connected to optical sources and 20 were connected to Si photodetectors [11,12].The fibers were placed at the mid-metatarsal level of the foot to provide images of a coronal cross section through the foot.This region was chosen to evaluate the vessels that supply the forefoot, the most common location for PAD foot ulcers.Our VOTI systems operated at two wavelengths (760 nm and 830 nm) and full tomographic data were obtained with a frame rate of approximately 3.3 Hz.
For the data acquisition, patients were in a sitting position with their foot placed inside the measuring system and a pressure cuff wrapped around their upper thigh.The sitting position was chosen instead of a standing position to avoid involuntary movements of the patients during the acquisition.A single measurement protocol consisted in 5 phases of exactly 1 minute each.It started with a resting phase (0:00 to 1:00 minute) in which a baseline signal was recorded; in the second phase (1:00 to 2:00 minutes) the cuff was inflated to 60 mmHg; in the third phase (2:00 to 3:00 minutes) the cuff was rapidly deflated and a new resting phase started to allow the signal to return to the baseline; in the fourth phase (3:00 to 4:00 minutes) a second cuff inflation was performed to 120 mmHg.Finally, in the fifth phase (4:00 to 5:00 minutes) the cuff was deflated, and the signal recorded till the end of the 5th minute.The aim during the cuff inflations was to obtain a venous occlusion so that the blood keeps being supplied to the foot through the arteries but its return to the heart is impeded.This results in a pooling of the blood in the foot, which in turn leads to increased light absorption.
The resulting data were input to a diffusion-theory-based PDE-constrained multispectral image reconstruction algorithm [27,28].For each patient, this algorithm calculated the total hemoglobin distribution for each frame acquired during the protocol, producing twodimensional cross-sectional images of the foot.We calculated the spatially dependent crosscorrelation coefficients ρ(h,W) between the area weighted average change in total hemoglobin over the whole cross section of the foot and the change in total hemoglobin for each pixel within the cross section over time: Here h is the time trace for an individual pixel in cross-sectional maps of the foot; W is the weighted average signal from the entire foot and µ h and µ W are their respective mean values in time.
A typical cross-correlation map ρ(h,W) for a PAD patient is shown in Fig. 3(A).It is possible to identify regions in which the local hemodynamic effects correlated strongly with the overall global hemodynamic effects.The "red" regions, in which the ρ(h,W) is high, correspond to the locations of the main vasculature in the foot.A comparison between the mean total hemoglobin variation in the entire cross section of the foot and just in the region were the correlation coefficient is higher than 0.7 (threshold value estimated in our previous study [12] that leads to the most significant differences between patients affected by PAD and healthy volunteers) is shown in Fig. 3(B).The changes in the regions of interest are considerably larger (up to ~22%) than changes observed when integrated over the entire foot (up to ~4%).Therefore, we used the same value here to obtain the experimental curves to compare with our model predictions.

Results and discussion
An example of a typical hemodynamic response during a cuff inflation/deflation (120 mmHg) is shown in Fig. 4 (black curve).Here the cuff is inflated at t = 10 seconds.Subsequently, we can see an increase in total hemoglobin for about 20 seconds, before a first plateau is reached.Sometimes, as in the case shown, this plateau becomes a "dip" in which the percentage variation of total hemoglobin (Δ[HbT]) briefly decreases, before a further increment is observed.A second maximum is reached at t = 70 seconds.At that point the cuff is deflated and Δ[HbT] returns to baseline.In addition to the measured hemodynamic response (black curve), Fig. 4 shows simulated data obtained from the physical model (green curves).Each green curve is obtained once the system (1) is solved as: represents the total blood volume (in milliliter) as the function of time in the area under observation by our optical imaging system.Both the experimental data and the simulated curves are then normalized using their maximum and baseline values.In this way we focus the attention on the behavior of the variation of milliliters of blood (and consequently total hemoglobin) during the cuff inflation/deflation, to reduce secondary effects depending on different foot/leg sizes that were not recorded during the previous experiment and could not be integrated in the model analysis.We found that certain physiological variables can be clearly connected to observed effects.For instance, the presence of the plateau or dip before the maximum value in the total hemoglobin variation (see Fig. 4) is strongly dependent on the R CUFF_Vm and especially its effect on the arteries expressed by the factor X.Only by varying these parameters can curves with two maxima and a dip in-between be generated.In Fig. 5(A), we varied X between 0 and 0.4 keeping all the other parameters to their minimum values reported in Table 1, and in Fig. 5(B) we varied R CUFF_Vm from 7 to 19 [mmHg s/mL] keeping X = 0.2.As can be seen, for X = 0 (which means the cuff does not affect resistance in the arteries) one observes only a steady rise of the signal towards a single maximum at 70 seconds.For X = 0.2 a maximum at 40 seconds starts to emerge, which becomes a pronounced peak for X = 0.4.Similarly, an initial maximum with varying peak values can be achieved by varying R CUFF_Vm .No other parameters in our model yielded similar effects.
Another important signal feature is the decline of the signal after 70 seconds (see Fig. 4).Here we found that the shape of this decline is strongly dependent on the resistances and capacitances of the veins (R V and C V ). Figure 5(C) shows results for 4 different values of R V ranging from 0.5 to 2 [mmHg s/mL] (always keeping all the other parameters to their minimum values reported in Table 1).Similar results can be obtained changing C V between 5 and 25 [mL/mmHg].Again, all other parameters had much weaker influence on that part of the signal.For example, Fig. 5(D) shows the effects of varying the resistance R A from 1.5 to 9.5 [mmHg s/mL].Only small shifts in the overall curve can be produced over the range of physiologically reasonable values (while keeping all the other parameters to their minimum values).After exploring the parameter space, we wanted to find the parameter that best represents the differences between healthy subjects and PAD patients.Therefore, we were looking for a parameter or sets of parameters that can be used to classify a hemodynamic response as belonging to a PAD patient with high sensitivity and specificity.Here we were guided by consideration of the underlying physiological changes known to occur in PAD patients.First, it is well established that plaque accumulation reduces the internal diameter of arteries [29] in PAD patients.According to Poiseuille law, the resistance is inversely proportional to the blood vessel radius to the fourth power [14,[18][19][20][21]. Consequently, one would expect a substantial higher value of R A and R AO in PAD patients compared to healthy subjects.
Second, it is known that deposit of lipids, cholesterol crystals, and calcium salt reduces the compliance of vessels in PAD patients.This means that the vessels do not expand in size as easily as vessels of healthy patient.Hence, during a venous occlusion, the vasculature of a PAD patients cannot accommodate as much blood as the vasculature of a healthy person [16,29]; or in other words, the capacitance of the vessel network is smaller.However, unlike the 4-th order effects of the reduction in vessel diameter, capacitance effects are only linearly proportional to the calcification, and its effects on the measurable signal is comparatively low [26,29].In addition, there is a significant variation in the vessel compliance due to aging [13,16,20,25], which is confounding the PAD diagnosis.
Based on these physiological considerations and our exploration of the potential parameter space performed above, we decided to use the total system resistance considering the cuff effect: as our biomedical marker to distinguish between patients with PAD and healthy volunteer.R T was determined for each of the subjects enrolled in our study by comparing the measured responses to the cuff experiment (see, e.g.Fig. 3) to 40 million curves in a look-up table we generated by varying the parameters in the range shown in Table 1.Therefore, for each experimental curve we found the set of resistance values that minimizes the square root of the absolute values of the differences between the experimental and model curve and we considered in the following analysis the average total resistance obtained from the best fittings of the two cuff inflations for each patient.
The mean values and standard deviation of R T were then calculated for the healthy volunteers and PAD patients.The results are shown Fig. 6.One can clearly see that the total resistance in healthy volunteers is considerably lower (27.1 ± 3.9 [mmHg s/mL]) than in PAD patients (33.9 ± 2.1 [mmHg s/mL]).An ANOVA analysis followed by a paired-Bernoulli test (using MATLAB embedded functions, anova1 and multcompare) shows that these differences are statistically significant between healthy and PAD patients (p value < 7 x 10 −6 ).Having established statistically significant differences in the R T values between healthy volunteers and PAD patients, we next performed a receiving operating characteristic (ROC) curve analysis (Fig. 7).Fig. 7. ROC curve to distinguish between healthy and PAD patients obtained considering the total system resistance when the cuff is inflated.
One can see that very higher sensitivities, specificities and accuracies are obtained when the total resistances are considered during the cuff inflation.For a cut-off value of R T = 31 [mmHg s/mL], we achieve a sensitivity of 95% and specificity of 95% (AUC = 0.97).
These results are however biased by the fact we used all our data to estimate the ROC curves, therefore, to further validate our findings, we performed a cross-validation analysis, based on a statistical resampling technique used for model assessment [21].To determine the optimal cut-off value for the total resistance to differentiate between healthy and PAD patients, a training subset of patients, randomly selected among all the patients, was used and then the remaining patients were used as testing samples to evaluate the model's generalization.This procedure was repeated changing the size of the training subset from 20 to 32 patients.Subsequently, the test subset changed from 20 to 8 patients, and, in each case, 500 reshuffling of the random combinations of patients to create the train and test sets were used to increase the objectivity of the assessment, obtaining the mean and standard deviations values showed in Fig. 8.In Fig. 8, we also applied the same cross-validation analysis to the percentage change in the total hemoglobin in the 60 mmHg peak only (which was used to obtain the results in [12]).It is possible to note that the estimated accuracy, sensitivity and specificity obtained with the model are higher than the ones obtainable by the percentage change in the total hemoglobin.We observed that, in the training, the mean values were always higher than 90% with standard deviations lower than 5% for all of them, whereas, in the test, accuracy and sensitivity had mean values around 90% with standard deviations lower than 10% and only the specificity had a standard deviation slightly higher (about 16%) and a mean value of about 86%.

Conclusions
In this paper, we introduced a new numerical model for simulating hemodynamic effects in feet, especially in response to a venous cuff occlusion.This model was tested on 20 healthy and 20 PAD patient measurements and provided a better understanding of the dynamic behavior of clinical data, obtained with vascular optical tomographic imaging systems, and lead to improvement of diagnostic sensitivity and specificity.
We found a correlation between the estimated total system resistance and the healthiness of the patient: the higher the total system resistance, the less healthy the patient.This was expected, since arteries in healthy subjects are less obstructed by calcifications, which reduce their internal diameter increasing the resistance to the blood flow, as compared to PAD patients.
By employing the new model, we determined mean values for sensitivity, specificity, and accuracy between 86 and 95% after a cross-validation analysis that showed a total standard deviation lower than 10% for accuracy and sensitivity and equal to 15% for specificity.These values compare favorably to other diagnostic technologies currently in use.For example, ABI can obtain high specificities (83.3-99.0%)and accuracies (72.1-89.2%),however at the cost of low sensitivities (15-79%) [6].Duplex ultrasound yields sensitivities between 80 and 85% and specificities higher than 95% [8], and studies on magnetic resonance imaging (MRI) have produced sensitivities between 85 and 90% and specificities of 95-97% [8].
In this paper, we limited the use of the model to estimate the variation of the total hemoglobin concentration to conform with the experimental data acquired.We did not test its ability to correctly estimate the absolute hemoglobin concentration in the patient.In future studies, in addition to the analysis of the blood accumulation in the tissues, information regarding the size of the subjects, their ages and weights will be recorded, and the model further validated on its ability to correctly estimate also the absolute hemoglobin concentration values.

Fig. 2 .
Fig. 2. Example of the time-dependent functions R V (t) (blue) and R A (t) (red) for a specific combination of the simulation parameters for a cuff inflation/deflation cycle.

Fig. 3 .
Fig. 3.A -The estimated cross-correlation coefficients ρ(h,W) in the cross sectional image of the foot calculated for [HbT] = [HbO2] + [Hb].B -The temporal map showing the percentage change of total hemoglobin distribution considering the entire cross-sectional area (red dotted line) and only the regions in which ρ≥0.7 (blue line).

Fig. 4 .
Fig. 4. Example of a typical hemodynamic response to a cuff inflation/ deflation (120 mmHg) recorded with our VOTI system (black curve).Also shown are 10 "best" fittings (green dotted lines) obtained automatically by the program.The cuff inflation starts at 10 seconds and lasts for 60 seconds.

Fig. 5 .
Fig. 5. Examples of the effects of X, R CUFF_Vm , R V and R A on the curves while keeping constant the other parameters.A-The effect of the cuff on the arteries (X) is the main cause in the creation of the experimentally observed plateau; B-Increasing the effect of the cuff (R CUFF_Vm ) creates a deeper plateau; C-The vein resistance (R V ) affects the curve only after the cuff starts to deflate; D-The main arteries resistance (R A ) effect is limited to a slight change in the slope when the other parameters are kept to their minimum values.

Fig. 6 .
Fig. 6.Mean and Standard Deviations of the average total system resistance during a cuff inflation estimated in Healthy and PAD patients.The mean values and standard deviations are respectively 27.1 ± 3.9 [mmHg s/mL] in Healthy and 33.9 ± 2.1 [mmHg s/mL] in PAD patients.

Fig. 8 .
Fig. 8. Train and Test results for both the model resistance (M) and the percentage change in total hemoglobin (PC) considering 500 different combinations of patients for each bar graph.The bar graphs are zoomed in the 60%-100% range to better visualize the changes in the mean values.Between 20 and 32 patients were used for the training (numbers on the left of the slash) and the remaining patients (from 20 to 8) were used for the test (numbers on the right of the slash).In black is shown the standard deviation obtained.