Methods for analyzing acetylcholine-induced vasodilation in human skin.

Data from healthy subjects from previously published (n=43 subjects) and unpublished (n=12 additional subjects) studies were included in the present analyses ( 29 , 34 – 39 ). The Institutional Review Board at The Pennsylvania State University approved all experimental procedures and protocols. Investigational Drug Numbers (IND 78,954; 124,294) were obtained from the Food and Drug Administration for the perfusion of N G -nitro-L-arginine methyl ester (L-NAME, non-selective NO synthase inhibitor; Calbiochem, EMD Millipore, Billerica, MA). Verbal and written informed consent were voluntarily obtained from all subjects before participation and in accordance with the guidelines set forth by the Declaration of Helsinki. Healthy men (n=9) and women (n=46) were screened for neurological, cardiovascular, metabolic, or dermatological diseases. The total data set included healthy control data from recent work in our laboratory, which has focused on disease specific to women (e.g., endometriosis and preeclampsia). Screening included a physical examination (heart rate, blood pressure, height, weight, etc.), medical health history questionnaire, and a blood chemistry analysis (Chem 24, Quest Diagnostics, Pittsburgh, PA). All subjects were not taking any prescription medications with known primary or secondary peripheral vascular effects. Women were tested without regard to the menstrual cycle phase or oral contraceptive use ( 40 ). The intradermal microdialysis technique has been described in detail ( 41 ). In brief, two microdialysis fibers (10 mm, 55 kDa, CMA Linear 31 probe, Harvard Apparatus, Holliston, MA) were inserted into the dermal layer of the ventral aspect of the forearm for the local delivery of acetylcholine (ACh, endothelium-dependent agonist; USP, Rockville, MD) dissolved in lactated Ringer’s alone or with 15 mM L-NAME (ACh+L-NAME). Perfusion rate was set to a constant 2 μL/min throughout the protocol (Bee Hive controller and Baby Bee microinfusion pumps; Bioanalytical Systems, West Lafayette, IN). After microdialysis fiber placement, a minimum of 60 min was allowed for resolution of the hyperemic response to needle insertion in all sites. During this period, L-NAME was perfused to allow for a drug wash-in phase. Duration of L-NAME perfusion is standardized, where prior work has demonstrated that a minimum of 60 min of L-NAME perfusion is sufficient to achieve NO synthase blockade and differences in baseline red blood cell perfusion between sites was not observed ( 20 , 42 ). Local heating units (VP12 and VHP2; Moor Instruments, Wilmington, DE) were set to a thermoneutral 33 °C directly over each microdialysis site ( 43 ). Laser-Doppler flowmeter probes were placed in the local heating units to measure local red blood cell flux (perfusion units). Following baseline measurements (~10 min), increasing concentrations of ACh from 10 −10 to 10 −1 M were sequentially perfused for 5 min at each dose to allow for steady state responses to occur in the ACh and ACh+L-NAME sites. Brachial blood pressure (Connex Spot Monitor, Welch Allyn, Skaneateles Falls, NY) was measured at the end of baseline and each ACh concentration perfusion period. After perfusing the final ACh dose (10 −1 M), maximal cutaneous vasodilation was induced by increasing local heater temperature to 43 °C and perfusing 28 mM sodium nitroprusside (SNP, NO donor; USP, Rockville, MD) in both sites ( 20 , 44 ). Blood pressure was obtained every 2–5 min until a plateau in red blood cell flux was obtained. Red blood cell flux data were recorded and stored for offline analysis (Windaq and Dataq Instruments, Akron, OH; PowerLab and LabChart, ADInstruments, Bella Vista, NSW, Australia). Average values for red blood cell flux were obtained during the following phases: minimum 5 min of baseline, minimum of the final 30 sec (and up to 2 min) of each ACh dose, and at the plateau during maximum vasodilation. Average red blood cell flux data were not obtained during the first minute of pharmacological agent perfusion to allow sufficient time for the agent to take effect. Mean arterial pressure (MAP) was calculated as the sum of one-third of the systolic blood pressure and two-thirds of the diastolic pressure at each phase. Absolute cutaneous vascular conductance (CVC) was calculated as shown in Equation 1 . (Eq. 1) CVC = Red Blood Cell Flux Mean Arterial Pressure To account for site-specific maximal vasodilation, relative CVC as a percentage of local maximal SNP and 43 °C heat ( % C V C m a x ) was calculated ( Equation 2 ) ( 45 , 46 ). (Eq. 2) % CVC max = C V C x C V C m a x × 100 C V C x was the averaged CVC measured at baseline or at a given ACh dose x and C V C m a x was the maximum averaged CVC for a given site. The robust regression and outlier removal (ROUT) method was applied to the absolute and relative CVC data to identify outliers independent of a model and its constraints, in which data were grouped by site and dose ( 47 ). The False Discovery Rate (Q) was set to 1% to identify outliers statistically while preserving biologically meaningful variability ( 47 ). Identified outliers were removed and the percentage of missing data was calculated as shown in Equation 3 . (Eq. 3) % Missing Data = # of Missing Observations Total # of Expected Observations × 100 Total number of expected observations is the product of number of subjects, sites, and doses including baseline and max. Outliers identified by the ROUT method, and missing observations, were removed and treated as missing values for all subsequent data and statistical analyses. Absolute CVC data were compared across site (ACh and ACh+L-NAME) and phase (baseline and max). A two-way repeated measures mixed-effects model was used with post hoc Tukey corrections (pharmacological site*phase, SAS 9.4, Cary, NC). Table results were presented as mean and standard deviation (SD). Statistical significance was set at α=0.05 in all instances. To assess pharmacodynamics, a four-parameter nonlinear regression model was used to generate group-averaged dose-response curves from individual %CVC max data for the ACh and ACh+L-NAME sites (Prism v10.5, GraphPad Software, La Jolla, CA). CVC max represents the maximal vasodilatory response to SNP perfusion and local heating at 43 °C. Therefore, CVC max within the relative CVC data ( Equation 2 ) were excluded. Best-fit parameters (top, bottom, logEC 50 , and Hill slope) with associated standard errors (SE) and their 95% confidence intervals (CI; based on the profile likelihood method) were derived for each curve ( 48 ). No constraints were applied and provided the following parameters: top, bottom, logEC 50 , and Hill slope. Poor model fit (i.e., failure to converge) indicates that the model’s predicted values did not fit the observed data and were identified in GraphPad by outputs including question marks, “unstable,” “very wide” parameter estimates, and/or extreme values out of the predefined axes limit (y = 100%; x = −1). Due to poor model fit, and thus high rate of subject exclusion, data analysis method 1 was further divided into four sub-analyses to assess sensitivity: 1A included all subjects regardless of best-fit values and 95% CI outputs, 1B included subjects with plausible (i.e., physiologically credible) best-fit values only, and 1C included all subjects with plausible best-fit values and 95% CI outputs. An ideal dose-response curve in an ex vivo experimental paradigm, which includes flat top and/or bottom plateaus, may not be consistently exhibited in in vivo models due to inherent variability. This violates the assumptions in nonlinear regression modeling, and thus, limits the ability of the data to be modeled. Therefore, 1D involved repeating and adding the value of the last ACh dose to the end of the dose-response curves to create an artificial top plateau. Subjects were included if plausible best-fit values and 95% CI outputs were generated. A fixed Hill slope of 1 (indicating no cooperativity) was implemented with no additional constraints. Data analysis method 2 provided the following parameters: top, bottom, and logEC 50 . Individual subject’s dose-response curves were excluded if the curves did not converge. Curves that did not converge were considered to have poor fit and could not be modeled using the nonlinear regression (curve fit) analysis. Top and bottom constraints with 100% and 0%, respectively, were implemented with no additional constraints. Data analysis method 3 provided the following parameters: logEC 50 and Hill slope. Individual subject’s dose-response curves were excluded if the curves did not converge. Following data analysis methods 1–3 and 1A-1D, the ACh and ACh+L-NAME dose-response curves and best-fit parameters generated within each method were compared using an extra sum-of-squares F-test (Prism GraphPad) ( 49 ). Methods 1–3 compared ACh and ACh+L-NAME curves and best-fit parameters within each method. For methods 1A-1D, sensitivity analyses were conducted to determine whether any approach altered the ACh and ACh+L-NAME curves or their associated best-fit parameters. F-test results are reported as the F-statistic with its associated degrees of freedom for the numerator and denominator, respectively [F (DFn, DFd)]. Table results are presented as best-fit values with associated SE and profile likelihood 95% CI. Figure results are presented as mean ± SD and 95% confidence bands. To assess site- and dose-dependent effects on microvascular endothelial function, group-averaged sites from individual %CVC max data at each dose within ACh and ACh+L-NAME sites were analyzed categorically using point-by-point analyses. CVC max data were excluded. A two-way repeated measures mixed-effects model was used to assess differences and post hoc Tukey corrections were used to account for multiple comparisons (pharmacological site*dose, SAS 9.4). Figure results are presented as mean ± SD. To quantify the total magnitude of microvascular response to the perfusion of ACh, the AUC of the %CVC max response to ACh was calculated using the trapezoid rule (Prism GraphPad) with the Y set to the averaged baseline. A baseline set to Y=0 is not physiologically plausible (implies zero skin blood flow) and artificially inflates the AUC value, which can overestimate the response that occurred following ACh perfusion. Moreover, these settings can capture the inherent variability of baseline CVC values among subjects, sites, and/or treatments. CVC max data were excluded. Minimum peak height and width were adjusted so that any peaks in the data were included in the analyses. Peak direction was also adjusted to include values that go above and below the baseline and total AUC was reported in arbitrary units (a.u.). A two-tailed, paired t -test was used to compare the AUC between ACh and ACh+L-NAME sites. The L-NAME-sensitive component of ACh-induced vasodilation, commonly referred to as NO contribution or NO-dependent vasodilation and a measure of microvascular function, was calculated as the difference between the AUC of the ACh and the ACh+L-NAME sites ( 24 , 34 , 50 ). Figure results presented as mean ± SD.

Subject characteristics are presented in Table 1 . Anthropometric, resting hemodynamics, and lipid profile characteristics were within normal clinical ranges. Absolute CVC data are shown in Table 2 . 3.9% of absolute CVC data was missing. Absolute baseline and max CVC were not different between sites (p=0.92). Absolute baseline CVC was significantly lower than absolute max CVC regardless of site (p < 0.0001). Relative CVC data were used for nonlinear regression (curve fit) modeling, point-by-point, AUC, and L-NAME-sensitive component analyses where the percentage of missing data was 2.6%. In data analysis method 1A, all subjects were included (ACh n=55, ACh+L-NAME n=55). ACh and ACh+L-NAME curves were significantly different [F (DFn, DFd): 133.6 (4, 1182), p < 0.0001; Table 3 , Figure 1A ]. LogEC 50 showed a rightward shift and Hill slope was reduced in ACh+L-NAME compared to ACh (p=0.0002 and p=0.02, respectively). There were no differences in best-fit values for top (p=0.36) and bottom (p=0.12) between curves. In data analysis method 1B, 12 and 14 subjects were excluded from the ACh and ACh+L-NAME curves, respectively (ACh n=43, ACh+L-NAME n=41). With the remaining subjects, ACh and ACh+L-NAME curves were different [109.1 (4, 899), p < 0.0001; Table 4 , Figure 1B ]. The top fit parameter was lower, and logEC 50 shifted rightward, in ACh+L-NAME compared to ACh (p=0.01 and p=0.0009, respectively). There were no differences in best-fit values for bottom (p=0.12) and Hill slope (p=0.06) between curves. In analysis method 1C, 26 and 38 subjects were excluded from the ACh and ACh+L-NAME curves, respectively (ACh n=29, ACh+L-NAME n=17). With the remaining subjects, ACh and ACh+L-NAME curves were different [71.78 (4, 497), p < 0.0001; Table 5 , Figure 1C ]. Best-fit value for logEC 50 showed a rightward shift in ACh+L-NAME compared to ACh (p=0.003). There were no differences in best-fit values for top (p=0.15), bottom (p=0.25), or Hill slope (p=0.16) between curves. In analysis method 1D, this approach yielded greater inclusion of data than 1C (ACh n=34, ACh+L-NAME n=24). With the remaining subjects, ACh and ACh+L-NAME curves were different [88.94 (4, 685), p < 0.0001; Table 6 , Figure 1D ; repeated final ACh dose plateau point corresponds to x=0]. The top parameter was lower and logEC 50 showed a rightward shift in ACh+L-NAME compared to ACh (p=0.0008 and p=0.0006, respectively). There were no differences in best-fit values for bottom (p=0.44) and Hill slope (p=0.27) between curves. Across approaches 1A-1D, ACh and ACh+L-NAME curves were not different [ACh: 1.52 (12, 1770), p=0.11; ACh+L-NAME: 1.54 (12, 1493), p=0.10]. Best-fit values for top (ACh: p=0.91; ACh+L-NAME: p=0.79), bottom (ACh: p=0.90; ACh+L-NAME: p=0.95), logEC 50 (ACh: p=0.12; ACh+L-NAME: p=0.32), and Hill slope (ACh: p=0.54; ACh+L-NAME: p=0.32) were not different within ACh and ACh+L-NAME curves across approaches, respectively. Out of 55 subjects, 1 and 5 subjects were excluded from the ACh and ACh+L-NAME curves, respectively (ACh n=54, ACh+L-NAME n=50). ACh and ACh+L-NAME curves were different [135.5 (3, 1121), p < 0.0001; Table 7 , Figure 2 ]. Best-fit values for top and bottom were lower in ACh+L-NAME compared to ACh (p < 0.0001 and p=0.02, respectively). There were no differences in best-fit values for logEC 50 (p=0.07) between curves. Out of 55 subjects, 4 and 12 subjects were excluded from the ACh and ACh+L-NAME curves, respectively (ACh n=51, ACh+L-NAME n=43). ACh and ACh+L-NAME curves were different [162.9 (2, 1018), p < 0.0001; Table 8 , Figure 3 ]. LogEC 50 shifted rightward and Hill slope was lower in ACh+L-NAME compared to ACh (both p < 0.0001). No subjects were excluded from this data analysis method (ACh n=55, ACh+L-NAME n=55). There were main effects of site and dose and an interaction effect (all p < 0.0001; Figure 4 ). Specifically, the %CVC max response at the ACh site was higher than ACh+L-NAME for a given ACh dose in the range of 10 −6 to 10 −1 M (all p < 0.01). No subjects were excluded from this data analysis method (ACh n=55, ACh+L-NAME n=55). The total AUC of the ACh site was greater than ACh+L-NAME [Mean (SD): 287.80 (89.29) vs. 167.40 (79.46) a.u., p < 0.0001; Figure 5A ]. The L-NAME-sensitive component showed a mean and SD of 120.30 and 112.70 a.u., respectively ( Figure 5B ). Individual values are shown in Figures 5A and 5B .

The utilization of the skin as an in vivo bioassay for assessing microcirculatory function has emerged as a unique integrative approach to examine physiological function in health and disease. This retrospective analysis presents five methods for analyzing cutaneous vascular conductance responses to a progressive ACh dose-response protocol (summarized in Figure 6 ). The primary findings of this study are that curve modeling methods for assessing ACh vasodilatory function in human skin yield different results depending on the method employed. Data analysis method 1 yielded the expected differences between curves, but had the highest rate of data exclusions and/or required more data reduction/processing (e.g., repeating an extra point at the top of the curve) of the data to fit the assumptions of traditional pharmacological curve modeling for inclusion of more data in the analysis. Data analysis methods 2 and 3 permitted the inclusion of more data with expected differences in the logEC 50 with data analysis method 3 allowing for detection of differences in the Hill slope. Using a mixed-effects model analysis (point-by-point) yielded a high degree of inclusion of data and provided insight into dose-specific changes between sites. Finally, analysis of the AUC and the L-NAME-sensitive component provided a more comprehensive picture into the overall vasodilatory response to increasing ACh doses, but with significant variability in the L-NAME-sensitive component. However, this analysis method provided additional insights into the underlying mechanisms involving NO synthase uncoupling. In data analysis methods 1–3, each curve modeling method confirms differences between the ACh and ACh+L-NAME curves and best-fit parameters (e.g., top, logEC 50 ). However, the number of exclusions of individual curves and the specific parameter differences were model-dependent. Data analysis method 1 exhibited the highest rate of exclusion. To explore this further, data analysis method 1 was analyzed in 4 different ways to limit excluded sites (data analysis methods 1A-1D) and tested for sensitivity. Dose-response curves and best-fit parameters (e.g., top, logEC 50 ) across data analysis methods 1A-1D were not different among the methods of data inclusion. However, specific differences in best-fit parameters between curves were model-dependent. Data analysis methods 1B and 1D were the only models able to detect the top parameter of the ACh+L-NAME curve was less than that of the ACh curve although all detected a rightward shift in logEC 50 . Given the large exclusion of data analysis methods 1C and 1D, methods 1A or 1B appear to be the most inclusive methods of analysis providing the most insight into the underlying physiology. For data analysis methods 2 and 3, as anticipated, best-fit values for top and bottom were lower in ACh+L-NAME compared to ACh. In data analysis method 3, best-fit values in logEC 50 were shifted rightward and the Hill slope was lower in ACh+L-NAME compared to ACh. If the intent is to explore changes in logEC 50 , then data analysis methods 2 and 3 are most appropriate with method 3 providing additional insights into the overall cooperativity of the response. We also analyzed the data set using a mixed-effects model (point-by-point) and AUC for ACh and ACh+L-NAME to quantify the L-NAME-sensitive component. Similar to data analysis methods 1–3, point-by-point analysis showed differences between ACh and ACh+L-NAME sites where %CVC max responses were higher in the ACh vs. the ACh+L-NAME site from ACh doses 10 −6 to 10 −1 M. Unlike data analysis methods 1B-1D and 2–3, point-by-point analysis allowed for inclusion of all subjects which reduced bias. In smaller data sets, using this analysis method would allow for the maintenance of statistical power. To account for the total response across doses to ACh, AUC analysis was used. As anticipated, the AUC for ACh+L-NAME was lower than the AUC for ACh. Interestingly, five subjects had a negative L-NAME-sensitive component (as low as −116.50 a.u.). In these subjects, the total vasodilatory response to ACh+L-NAME was increased over that of the ACh site. While more commonly observed in subjects with overt endothelial dysfunction, a negative L-NAME-sensitive component can provide insight into how endothelium-dependent mechanisms may be altered. A plausible physiological explanation for a negative L-NAME-sensitive component includes the inhibition of the production of oxidant species through uncoupled NO synthase. These specific mechanisms require further interrogation in both clinical and healthy populations. Endothelial dysfunction measured in the cutaneous vasculature evolves in parallel and may precede overt detectable changes in vascular dysfunction to other circulatory beds ( 5 – 10 ). Debbabi and colleagues reported a close relation between microvascular endothelial function and conduit artery endothelial function following ACh-mediated vasodilation ( 51 ). Their results show the integrated contribution of several mechanisms such as NO, COX, and EDHF to the endothelial response to ACh. In conjunction with ACh+L-NAME, where NO synthase activity is non-selectively inhibited, the role of NO relative to other endothelium-dependent pathways can be further delineated. We and others have used traditional pharmacological curve modeling analyses in our in vivo skin bioassay. There are several advantages to this type of analysis from a physiological context including being able to discern differences in the sensitivity of the cutaneous vasculature to specific endothelium-dependent agonist. Additionally, this approach aligns closely with preclinical pharmacologic methods, highlighting its translational utility. However, in practice, pharmacological curve modeling in the in vivo model often does not result in the best-fit of the individual data and key assumptions are not met. There are several examples of the assumptions for pharmacological curve modeling being violated using this in vivo model. ACh perfusion through intradermal microdialysis in healthy vasculature can occasionally induce higher flux values at the 10 −2 and 10 −1 concentrations in comparison to the flux response induced with high dose SNP (28–50mM) and local heating at 43 °C. In these cases, it is best practice to use the highest flux (CVC) value as the CVC max when converting absolute CVC to relative CVC ( Eq. 2 ). This response likely reflects a degree of tachyphylaxis from the vessels being exposed to excessive vasodilatory stimuli. Further, in ex vivo models, baseline is standardized using pre-constriction ( 52 , 53 ). Expressing CVC as relative to site-specific maximum, baseline standardization considerations, and top and bottom constraints (data analysis method 3) produces best-fit estimates for logEC 50 (sensitivity) and Hill slope (cooperativity). A Hill slope constraint of 1 (data analysis method 2) allows the assumption of a standard, non-cooperative dose-response relation where the vasodilatory response to ACh is concentration-dependent as opposed to reliant on receptor cooperation. This constraint produces best-fit estimates for logEC 50 and minimal/maximal vasodilatory responses. Taken together, the implementation of curve modeling methods should be dependent on the underlying physiological question. Our data suggests that the unconstrained model (data analysis method 1A or 1B) provided curves and best-fit values that most accurately reflected the observed responses at each dose and site, whereas the constrained fits deviated substantially. Due to the variability in an in vivo model, and the necessity for site exclusion with curve modeling, the point-by-point analysis, with AUC and L-NAME-sensitive component calculations, quantifies the contribution of NO and limits the potential exclusion of data points. The majority of the subjects used in this retrospective analysis were female. Menstrual cycle and oral contraceptive use were not controlled ( 40 ). The impact of menstrual cycle and hormonal contraceptive methods on microvascular function has been the topic of considerable debate, but data remain equivocal in healthy women ( 54 ). Here, we offer best practices for analyzing and reporting cutaneous vascular conductance responses to a progressive ACh protocol. First, utilize the analytical approach that most accurately addresses the physiological question. If the purpose of an investigation is to assess and quantify NO-dependent vasodilation, point-by-point, AUC, and L-NAME-sensitive component analyses are recommended, especially when dealing with smaller sample sizes. The use of nonlinear regression (curve fit) modeling can detect whether differences exist between curves, as well as curve parameters to provide pharmacodynamic information (e.g., logEC 50 , Hill slope), but this approach is limited by data inclusion and cannot directly assess NO-dependent vasodilation. If nonlinear regression (curve fit) modeling is utilized, data analysis methods 1A and 1B are recommended as the curves both reflect the visualized data with comparable accuracy to methods 1C and 1D. Methods – including how outliers are defined, how missing data are handled (if at all), models and constraints used, and what values (i.e., group-averaged from individual %CVC max data) are reported – should be clearly stated to strengthen rigor and reproducibility. Accordingly, all parameters generated by the models and constraints should be reported in full rather than limiting the report to the parameter(s) of interest. Reporting accurate goodness-of-fit metrics with nonlinear regression models should always be done. For instance, reporting of R 2 as it relates to nonlinear regression modeling should be avoided as it can be misleading – R 2 assumptions are more appropriate for linear regression models ( 55 ). Furthermore, goodness-of-fit metrics such as standard error and 95% confidence intervals should be reported as both aid in evaluating and providing insight into parameter precision and model stability. In addition, AUC calculations should be calculated as described above and detailed to ensure that investigators can make group and/or treatment comparisons across various labs. This alone can avoid masking or exaggerating treatment effectiveness and therefore improve clarity and accurate interpretation of the results. Earlier work has discussed the use of conductance versus resistance in blood flow related work ( 56 – 58 ). Conductance (flow/pressure gradient) is the inverse of resistance (pressure gradient/flow) and both are used as indices of vascular tone ( 56 ). Though the relation to blood flow when the pressure gradient does not change is different: conductance has a linear, and resistance a nonlinear (hyperbolic), relation to blood flow. Reporting of conductance is the preferred metric, as any change to blood flow is proportional to conductance. Conversely, since resistance has a hyperbolic relation with blood flow, the magnitude of change is dependent on where one is on the hyperbolic curve – either significant or minimal changes can occur which can distort the findings and interpretation of the data ( 56 – 58 ). In the context of aging or pathophysiology, careful consideration for reporting absolute and/or relative CVC should be considered. As demonstrated in the present paper, absolute baseline and max CVC were compared between sites prior to the use of relative CVC. In the event where differences in phases (baseline or max) are observed across sites, treatments, or groups, absolute CVC should be used to make comparisons, even though absolute CVC may introduce greater variability due to site-specific differences in the number of microvessels beneath the laser-Doppler probe. Further, when key group differences exist by design such as in comparing normotensive and hypertensive samples, it may be necessary to report both raw flux (i.e., index of flow) in addition to conductance because there are differences in overall driving pressure ( 24 ). We suggest reporting data from the progressive ACh dose protocol to use EC 50 as opposed to ED 50 or IC 50 . EC 50 is the concentration of the agent that produces the half-maximal response, whereas ED 50 is the concentration of the agent that produces a specific effect in 50% of the population that has been administered that dose ( 59 , 60 ). IC 50 , on the other hand, is the concentration of the agent required to inhibit a response by half ( 25 ). All in all, EC 50 is the term that is most applicable to the data. ACh concentrations are also commonly log-transformed, so reporting logEC 50 is most accurate. It has also been common practice to designate the difference between the ACh and ACh+L-NAME AUCs as NO contribution or NO-dependent vasodilation; it is more precise to refer to this as the L-NAME-sensitive component. Under pathophysiological conditions, blockade of NO synthase with L-NAME may alter oxidant generation through uncoupled NO synthase and thereby increase the overall vasodilatory response to ACh. This would result in a negative L-NAME sensitive component to ACh-mediated vasodilation. These negative values are often unreported as they may be viewed as spurious data points. Therefore, it is highly encouraged to report such findings with individual data points to improve transparency.

The goal of this retrospective analysis was to report on five methods for analyzing cutaneous vascular conductance responses to a progressive ACh protocol. All nonlinear regression (curve fit) modeling methods showed differences between ACh and ACh+L-NAME curves. Data analysis methods 1A and 1B yielded curves and best-fit values most consistent with visualized data. Although nonlinear regression modeling can provide pharmacodynamics of ACh and ACh with L-NAME in an in vivo human model, the analysis methods may pose limitations if the data cannot be modeled due to poor fit. Point-by-point analysis, and calculations of the AUC and the L-NAME-sensitive component, allow for direct quantification of NO-dependent vasodilation, but are unable to characterize pharmacodynamic responses. However, the latter analysis methods circumvent the necessity for data exclusion often observed with curve modeling and may be more appropriate for small sample sizes.

Endothelial dysfunction is characterized as impaired endothelium-dependent vasodilation, augmented vasoconstrictor responsiveness, and vascular remodeling ( 1 , 2 ). Reduced bioavailability of endothelial-derived vasodilators, including nitric oxide (NO), are implicated in the development of cardiovascular disease ( 3 , 4 ). Endothelial dysfunction occurs in parallel across different vascular beds (e.g., renal, coronary, and skin circulations) appearing early in the microcirculation and predicting long term adverse cardiovascular outcomes ( 2 , 5 – 12 ). The human cutaneous circulation is an accessible vascular bed used to directly assess mechanisms of microvascular function and dysfunction ( 1 , 8 ). The dense capillary network in the skin is ideal for delivery of pharmacological agents to investigate neurovascular signaling mechanisms ( 1 , 13 – 16 ). Intradermal microdialysis is a minimally invasive technique that allows bidirectional delivery of small molecular weight substances directly to the cutaneous microvasculature ( 13 , 14 ). Integrating intradermal microdialysis with laser-Doppler flowmetry and local heating allows for measurement of red blood cell flux changes (i.e., vasoreactivity) to localized perfusion of pharmacological agents. The localized nature of drug delivery with intradermal microdialysis prevents systemic effects and provides a level of within-person control that is advantageous over other techniques including intravenous or intra-arterial infusion ( 17 ). Acetylcholine (ACh) is a receptor-mediated, endothelium-dependent agonist that binds to muscarinic receptors and induces vasodilation through various second messenger pathways including NO, cyclooxygenase (COX), and other endothelium-derived hyperpolarizing factors (EDHF) ( 18 – 20 ). In the skin, ACh dose-response protocols in which ACh is dissolved in a vehicle (e.g., saline, lactated Ringer’s) and perfused alone (ACh) or concurrently with N G -nitro-L-arginine methyl ester (ACh+L-NAME, a non-selective NO synthase inhibitor) are used to directly measure endothelium-dependent NO-mediated vasodilation. This in vivo bioassay has been extensively used to examine physiological and pathophysiological mechanisms mediating vascular dysfunction in the context of aging, hypertension, diabetes, and depression among others ( 19 – 24 ). Although dose-response curves are commonly used to describe in vivo cutaneous vasodilator responses, modeling parameters and analytical methods are not consistent across studies. A common analysis method is nonlinear regression (curve fit) modeling, which provides pharmacodynamic information ( 25 , 26 ). Curve parameters obtained from the nonlinear regression include best-fit parameters such as: top (maximum vasodilator response), bottom (minimum vasodilator response), logEC 50 (concentration of the agonist, on a log scale, that elicits 50% of the maximal response), and Hill slope (degree of cooperativity between the agonist and receptor) ( 25 ). However, best-fit parameters generated from group averaged curves, or averaged estimates from individual fits, and constraints placed (e.g., variable or fixed Hill slope) are inconsistently reported ( 22 , 27 – 30 ). Alternatively, point-by-point analysis involves statistical comparisons of vascular responses across ACh concentrations (i.e., doses) and between different local treatments ( 19 , 22 , 27 – 29 , 31 – 34 ). Calculating the differences between area under the dose-response curves (AUC) for ACh and ACh+L-NAME sites allow for quantification of the total magnitude of microvascular response – the L-NAME-sensitive component to ACh-mediated vasodilation ( 27 – 29 , 34 ). In both curve-modeling and point-by-point analyses, discrepancies remain in how missing data and AUC calculations are performed ( 19 , 27 , 29 , 32 , 33 ). To date, there are no established recommendations for analysis of skin microcirculatory responses to vasodilators. Therefore, the aim of this retrospective analysis was to evaluate different methods for analyzing microcirculatory responses to ACh in human skin. These methods included nonlinear regression (curve fit) modeling with 1) no constraints, 2) Hill slope constraint, and 3) top and bottom constraints; as well as 4) point-by-point analysis, and 5) AUC calculation for ACh and ACh+L-NAME sites with the calculated L-NAME-sensitive component.