Heterogeneity in patterns of pain development after nerve injury in rats and the influence of sex.

The genesis of neuropathic pain is complex, as sensory abnormalities may differ between patients with different or similar etiologies, suggesting mechanistic heterogeneity, a concept that is largely unexplored. Yet, data are usually grouped for analysis based on the assumption that they share the same underlying pathogenesis. Sex is a factor that may contribute to differences in pain responses. Neuropathic pain is more prevalent in female patients, but pre-clinical studies that can examine pain development in a controlled environment have typically failed to include female subjects. This study explored patterns of development of hyperalgesia-like behavior (HLB) induced by noxious mechanical stimulation in a neuropathic pain model (spared nerve injury, SNI) in both male and female rats, and autonomic dysfunction that is associated with chronic pain. HLB was analyzed across time, using both discrete mixture modeling and rules-based longitudinal clustering. Both methods identified similar groupings of hyperalgesia trajectories after SNI that were not evident when data were combined into groups by sex only. Within the same hyperalgesia development group, mixed models showed that development of HLB in females was delayed relative to males and reached a magnitude similar to or higher than males. The data also indicate that sympathetic tone (as indicated by heart rate variability) drops below pre-SNI level before or at the onset of development of HLB. This study classifies heterogeneity in individual development of HLB and identifies sexual dimorphism in the time course of development of neuropathic pain after nerve injury. Future studies addressing mechanisms underlying these differences could facilitate appropriate pain treatments.

Introduction

While all neuropathic pain disorders are the result of damage of the somatosensory system, their etiologies are diverse, ranging from trauma and diabetes to genetic disorders, exposure to cancer chemotherapy, or infections. (Costigan et al., 2009, vonHehn et al., 2012). The complex genesis of neuropathic pain is also evident from observations that sensory abnormalities also may differ between patients with similar neuropathic etiology, suggesting additional unexplored mechanistic heterogeneity (Baron et al., 2017, Baron et al., 2010, Maier et al., 2010). Variability is likewise evident between individual animals in preclinical models of neuropathic pain (Dean et al., 2017, Hogan et al., 2004, Jaggi et al., 2011, Roytta et al., 1999). Despite these observations, data are usually grouped for analysis based on the assumption that they share the same underlying pathogenesis. An alternative approach is to consider that animal subjects may be members of distinct pathogenic groups despite their apparent baseline similarities. In this study, we examine animals subjected to identical nerve injuries but seek to identify distinct groups based on each animal’s time course for the development of hyperalgesia.

Sex is a factor that may contribute to individual differences in pain responses, as the incidence of chronic pain is considered greater in women (Berkley, 1997, Gerdle et al., 2008, Unruh, 1996, Wijnhoven et al., 2007), who generally have higher sensitivity and a lower tolerance to pain (Dixon et al., 2016, Mogil, 2012, Riley et al., 1998), and report higher pain ratings than men (Ruau et al., 2012). However, these conclusions are inconsistent, often drawn from underpowered studies (Berkley, 1997, Mogil, 2012, Racine et al., 2012, Riley et al., 1998). Furthermore, sex differences are typically variable across pain modalities and measures, and often conclusions are biased by highlighting a single statistically significant measure of a given outcome to the exclusion of others that were not (Mogil, 2012, Racine et al., 2012). Trauma of the nervous system often results in neuropathic pain presenting as hyperalgesia and/or allodynia (Woolf, 1995) and in the clinical realm, is considered more prevalent in women (Bouhassira et al., 2008, Colloca et al., 2017, de Mos et al., 2007, Sandroni et al., 2003, Torrance et al., 2006). However, female subjects are inadequately studied in preclinical research (Mogil, 2012, Mogil and Bailey, 2010), and the limited findings are conflicting. Female rats have shown greater susceptibility to developing neuropathic pain (Coyle et al., 1995) and more intense allodynia, but this is dependent upon strain (DeLeo and Rutkowski, 2000, LaCroix-Fralish et al., 2005, Nicotra et al., 2014). Other studies in rats have found no sex differences in withdrawal threshold or mechanical hypersensitivity after sciatic nerve injury (Dominguez et al., 2009), chemotherapy-induced neuropathy (Legakis et al., 2018), partial spinal nerve ligation (Severino et al., 2018) or plantar foot incision (Kroin et al., 2003), with similar findings in a mouse model for chemotherapy- and SNI-induced mechanical allodynia (Bourquin et al., 2006, Naji-Esfahani et al., 2016). While sex differences in mechanical and cold allodynia have been assessed in the SNI model, no study has reported sex differences in pin-prick hyperalgesia, a test stimulus that is aversive, unlike the von Frey test (Wu et al., 2010). In addition, much consideration has been given to quantitative differences in pain intensity, but less to temporal differences where pain end-point is the same, but is driven by different mechanisms (Juni et al., 2010, Mogil, 2012, Sorge et al., 2011, Tershner et al., 2000), a concept that may apply to neuropathic pain (Severino et al., 2018). For these reasons, we have included both sexes in this examination of underlying diversity in the development of pain.

It is well established that the somatic sensory system and cardiovascular function are closely linked to maintain homeostasis, (Keay and Bandler, 2001) and that sympathetic function is altered in some patients with chronic pain (Bruehl and Chung, 2004, Nordin and Fagius, 1995, Schlereth and Birklein, 2008, Walters, 2018). Our prior study in male rats demonstrated that development of hyperalgesia-like behavior (HLB) after nerve injury is associated with a decrease in sympathetic tone (Dean et al., 2017). These observations suggest that maladaptive changes in sympatho-sensory regulation could contribute to chronic pain conditions and autonomic disturbances (Lovick, 1990, Millan, 2002). While the primary goal of this study was to extend our prior work (Dean et al., 2017) to explore heterogeneity in development of neuropathic pain and assess the influence of sex in gonadally intact animals in a manner that considers the full time course of post-injury changes, a secondary objective was to explore the interaction of pain behavior with autonomic function to provide insight into the integrative nature of pain. For these studies, a well-documented spared nerve injury (SNI) procedure which produces hypersensitivity to mechanical stimuli in the rat was used (Dean et al., 2017, Decosterd and Woolf, 2000).

This analysis used several approaches to categorize the time course of change in HLB and autonomic function after nerve injury, and the influence of sex. Specifically, data were longitudinally clustered using discrete mixture modeling and using a new rules-based pain trajectory classification scheme. In analysis, mixed models were used to identify linear relationships because they accommodate correlations and differing variances within the data, and they have been shown via simulation to detect differences with smaller sample sizes with the same level of power than other methods (Aschenbrenner et al., 2014).

Materials and methods

The protocols for the study were approved by the Animal Care and Use Committees at the Medical College of Wisconsin and the Zablocki Department of Veterans Affairs Medical Center. Male (n = 82, 255–375 gm) and female (n = 79, 249–348 gm) Sprague-Dawley rats were obtained from Charles River (Wilmington MA) and were maintained and used according to the NIH Guide for the Care and Use of Laboratory Animals, and in compliance with federal, state, and local laws. This study is in what is considered the first stage of testing for a sex difference, which examines gonadally intact adult male and female rats (Greenspan et al., 2007, Prendergast et al., 2014). Addressing the influence of the estrous cycle is beyond the scope of this paper and testing the specific role of sex hormones will be the focus of further study. Although single housing could influence social aspects of pain, all on-study animals were housed individually to avoid confounding, as individual housing is necessary for radiotelemetric blood pressure monitoring. Animals were housed and studied undisturbed in a room on a reverse 12 hr light/12 hr dark cycle maintained at 22 ± 1 °C at 35–45% humidity. Animals had free access to food (Purina laboratory rodent diet 5001) and water, and bedding was Beta Chip (Warrensburg NY).

Instrumentation and blood pressure monitoring

A subset of animals (n = 45 male; n = 44 female) were implanted with PA-S10 or PA-X10 transmitters (Data Sciences International (DSI), St. Paul, MN) for radiotelemetric monitoring of arterial blood pressure. Telemetry surgery was performed 7 days prior to SNI surgery. Rats were anesthetized with isoflurane (5% induction, 1.75–2% maintenance) in O2, and the cannula (0.43 mm outer diameter polyethylene) attached to the transmitter was inserted into the left femoral artery through a small inguinal incision, with care taken to avoid manipulation of the adjacent femoral nerve. The transmitter was fixed in a subcutaneous pocket on the left flank of the rat and the incision was closed with 3–0 silk suture. A redundant loop in the cannula allowed for the growth of the rat. Animals were treated with carprofen (5 mg/kg, sc) at the start of surgery and 24 h postoperatively. After 7 days, animals were again anesthetized for removal of inguinal sutures and for SNI or sham surgery (details below).

Spared nerve injury (SNI)

All spared nerve injuries were performed by the same experienced surgeon who has performed this surgery regularly (Dean et al., 2017). Training for this surgeon included development of reliable identification of individual nerves using electrical stimulation to trigger expected motor responses. An incision was made through the skin and the biceps femoris muscle on the lateral surface of the right thigh to expose the sciatic nerve and its three terminal branches (tibial, common peroneal, and sural nerves). The tibial and common peroneal nerves were ligated with 6.0 silk suture and transected, while care was taken to avoid trauma to the sural nerve, which remained intact. The muscle layers were closed with 5-0 absorbable polyglycolic suture and the skin was closed with staples that were removed after 7 days. Carprofen (5 mg/kg sc) was administered at the start of surgery. No postoperative analgesia was administered to avoid confounding effects on hypersensitivity or development of pain behaviors post injury.

Sham SNI surgery was performed on an additional control group of female rats (n = 4), in which all three terminal nerves were exposed but left intact. Experimenter was not blinded to this being a control group. We have previously published data for sham control male rats that demonstrates an absence of HLB (Dean et al., 2017).

Data analysis

Behavioral evaluation of HLB

Mechanical hypersensitivity, a measure of HLB, which is an enhanced response to a noxious stimulus was evaluated as previously described after each resting blood pressure monitoring period (Dean et al., 2017, Rigaud et al., 2011). Animals were placed individually in clear plastic enclosures on an elevated ¼ in wire grid for a 15 min acclimation period that allowed animals to cease exploratory activity. The point of a 22 g spinal anesthesia needle was applied to the lateral part of the glabrous plantar surface of the paw with sufficient needle force to indent, but not penetrate, the skin. The behaviors induced by this noxious stimulus were of two types, either a brisk, simple withdrawal with immediate return of the foot the wire floor, which is typical of normal animals, or a hyperalgesia-type response that consisted of sustained elevation of the paw with shaking, licking and grooming, for which we use the term HLB (Gemes et al., 2009, Hogan et al., 2004). Although there are numerous measures of an animal’s pain experience, this test was chosen because hyperalgesia induced by a pin touch is a common finding in clinical pain patients, especially those with peripheral neuropathic pain (Bennett, 2001, Rasmussen et al., 2004, Scholz et al., 2009). Additionally, we have confirmed that positive responses to this test specifically correlate with conditioned place avoidance in plantar skin testing of nerve-injured rats, thus identifying an aversive experience, whereas other tests such as simple response to calibrated monofilaments to determine withdrawal threshold (von Frey test) are not aversive (Wu et al., 2010). Conditioned place avoidance cannot be repeated, and as such is not suited for a longitudinal study design as in this project. Therefore, the frequency of HLB responses during needle testing was used as the most meaningful measure of pain hypersensitivity. One individual assessed the response type for each of 5 applications to each hindpaw, ipsilateral and contralateral to the injury. Mechanical stimuli were separated by at least 10 s, and repeated after 5 min, for a total of 10 touches to each paw. For each time point, the total number of HLB responses (a maximum of 10) was converted to % to provide the HLB response rate for analysis. In agreement with our previous findings, all rats responded with a brisk withdrawal of the paw but no HLB on stimulation of the contralateral paw, all reported data were on the ipsilateral paw (Dean et al., 2017). Baseline HLB response rate prior to SNI was 0% for all but 8 (5.0%) animals and was not different by sex (p = 0.2865). In addition, baseline HLB response rate for animals with telemetry surgery was the same as for those without (p = 0.0988). All animals underwent SNI except the separate group of 4 Sham control females, and therefore condition was unblinded, as was sex. The same investigator evaluated all behavior and was blind to the autonomic findings of the subjects.

Heart rate variability determination

In all animals, blood pressure was monitored at least one hour prior to behavioral evaluation of hyperalgesia to avoid confounding effects of sensory testing. No one was present for monitoring which took place in the room in which the animals were housed so that they did not need to be moved or require a period of habituation. Resting blood pressure was analyzed for heart rate and heart rate variability (HRV) to indicate autonomic changes to the heart, using the DSI Ponemah 6.4.2 acquisition, analysis and archive system. Spectral analysis of heart rate was used to determine HRV from blood pressure data as described previously (Gemes et al., 2009). The blood pressure data were used to calculate an inter-beat interval (IBI) series from three 5-min segments, each with 5 overlapping (50%) sub-segments of 512 points. The IBI was converted into an instantaneous heart rate using the formula: heart rate (beats/min) = 60/IBI (sec). The IBI values were then interpolated at 50 Hz (cubic algorithm), detrended, and the mean was suppressed to create the data points equally spaced in time for further analysis. Utilizing a Hanning window, the power was calculated for each data set over the frequency ranges of 0.25–1 Hz (low frequency, LF) and 1–3 Hz (high frequency, HF), with the results from the three segments averaged for the final determination of density within each frequency band. LF and HF powers were normalized and expressed as percentages of total power (LF + HF powers) and the LF/HF ratio was calculated from these values as a representation of sympathetic tone from sympatho-vagal balance of direct autonomic influence on the heart (Eckberg, 1997, Malliani et al., 1991).

Behavioral test/blood pressure monitoring schedule

HLB response rate and LF/HF power ratio of HRV were calculated for each monitoring period on days 1, 3, 7, 10, 14, 18 and 21 post injury. Baseline (pre-nerve injury, day -1) levels were calculated by averaging values for the 2 days prior, and the day of but prior to, nerve injury. The 4 sham female rats were tested and monitored on the same schedule.

Analysis of patterns of behavior

The analysis plan was to initially explore the raw data graphically, blinded to sex. These data indicated that temporal pattern assessment was necessary because successive pain measurements are expected to be related. Consequently, two methods of temporal classification with different strengths were used and comparisons across sex were made after classification.

Temporal patterns of behavioral phenotype were identified according to the sequence of changes in hyperalgesia response rate over the 21 days post-injury. Baseline responses were not used because the interest was in the development of post-SNI behavior. A Moving Average filter with lag-1 (MA(1)) was applied to reduce random volatility of the hyperalgesia response data, by averaging each two consecutive data points, this was not done with the sympathetic data because the volatility (variation) is of interest by definition. Missing values in both data types were replaced by the average of their preceding and following data points (HLB response data were not collected for 23 males, and HRV not assessed for 22 of those males, on day 18. BP data could not be analyzed for HRV in 6 males and 1 female on day 1, 1 male and 1 female on day 7, 2 females on day 10, and 6 females on day 18), rats with two consecutive missing values were excluded from analysis (hyperalgesia response data was not collected for 14 males on days 1 and 3. BP could not be analyzed for HRV in 1 female on days 10 and 14, 1 female on days 14 and 18, and 3 females on days 18 and 21). Time course patterns were defined using two different methodologies for the hyperalgesia data, seeking concurrence in conclusion. The longitudinal discrete mixture modeling method evaluated similar temporal sequences, primarily identifying groups according to data magnitude. In contrast, the rules-based method was defined a priori, and identified groups by overall trend, scaled by magnitude.

While the results of mixed models and repeated measures ANOVA have different interpretations, both identify whether temporal relationships are present (Aschenbrenner et al., 2014). Because the latent-class groupings were defined by magnitude, mixed models were used to fit their relationship with time, because the rules-based groupings were defined effectively a priori, repeated measures ANOVA was used to determine whether there was a temporal relationship.

Discrete mixture model longitudinal clustering.

The discrete mixture model clustering method, applied to both HLB and HRV outcomes, sought groups of cases with similar values in the same temporal sequence, and so primarily identified groups according to data magnitude. Thus, for each group, plotted y-axis values show the expected magnitude, based on the contributing data points, over time. Rats were grouped via TRAJ, a SAS macro (SAS Institute, NC) that performs longitudinal clustering semi-parametrically using discrete mixture modeling (Jones et al., 2001). Because this is a method of latent class identification, it is less influenced by the analysts preconceived notions, although it is acknowledged that to obtain a result, with high positive predictive value, a sample of >450 is required, which is vastly different from sample sizes of <20 typical used in preclinical research.

Clustering was performed with both sexes combined, and then separately for each sex. Unless stated otherwise, the reported results are from the clusters derived from the combined sex data. The ideal number of clusters for each dataset was identified using 2(ΔBIC), as recommended (Jones et al., 2001).

Rules-based longitudinal grouping.

Rules-based longitudinal grouping was defined a priori based on knowledge of pain development and mathematically distinct line shapes and identified groups of cases by shape (overall trend), scaled by magnitude. This method was only applied to the HLB data. Rats were grouped according to a set of inclusive, explicit, and non-overlapping categories established independently from the data and without influence of sample size. Changes to overall shape were measured via differences between temporally consecutive points; they were considered different if the difference was ≥25% of one of the adjacent data point/s (as specified below). This value was set prior to examining the data, to exceed a HLB response rate of 20% which is selected to be outside of the range of random error (Dean et al., 2017, Gemes et al., 2009). The rules were defined after data collection.

Nine behavior categories were defined as below and shown in Table 1, and for secondary analysis, the two increasing patterns (BH) were combined, as were the two decreasing patterns (CG). Assignment of category was performed twice, with a time gap of more than one month to avoid recall bias, by an analyst blinded to the sex of the animal and its other outcomes (kappa = 0.9384).

Table 1

Rules-based longitudinal grouping of hyperalgesia behavior after SNI by sex.

Zero, All values = 0

Non-decreasing, allowing for one negative difference of < 25% of the temporally earlier data point

Non-increasing, allowing for one positive difference of < 25% of the temporally later data point

Hump-shaped – strictly non-decreasing, then a maximum (more than one consecutive time-point may hold the maximum value), then strictly non-increasing (the differences between either of the maximum’s two immediately preceding data points and either of the maximum’s two immediately following data points must both be ≥ 25% of the maximum)

U-shaped – same parameters as D but in the opposite direction

Sinus-shaped – displays both a hump- and a U- shape such that the maximum of the hump and the minimum of the U are not immediate data points; the hump and U may appear in either order

For cases not satisfying the criteria of groups A-F, the first two and last two values in the temporal sequence were examined. If the difference between each of the pairs of values was ≥ 25% of the smaller (but non-zero) value present, the more extreme value (the value further from the values of the opposite pair) was retained as the representative. Otherwise, the more moderate value (the value closer to the values of the opposite pair) was retained as a representative for further assessment.

Otherwise decreasing – the difference between the MA(1) starting and finishing representatives was ≥25% of the larger value and the starting representative was >the finishing representative

Otherwise increasing – the difference between the MA(1) starting and finishing representatives was ≥25% of the larger value and the starting representative was

Otherwise flat – the difference between the MA(1) starting and finishing representatives was <25% of the larger value

Statistical analysis

The study was originally powered for two groups, those that developed a HLB response rate of 20% on day 21 and those that did not. Review of this initial grouping (40 males, 50 females) showed a sex difference signal that could not be confidently interpreted because the statistical methods (ANOVAs) were too coarse, which in turn rendered the sample size too small. Additionally, the temporal patterns prior to day 21 indicated that a binary grouping lacked sensitivity to the primary concern of HLB-development. To increase sample size but reduce animal use, some data for the 21 days post-injury for this report was acquired from protocols that pursued other avenues of investigation after the 21 day monitoring period reported here. All animal care, monitoring and testing protocols for the pre and 21-day post injury periods were identical. Consistency was achieved with SNI performed by the same surgeon and all behavior tests carried out by the same individual. The goal was to obtain similar numbers of males and females.

Analysis was performed using SAS 9.4–14.3. Visualizations were performed in R 3.6.0. Data grouped via discrete mixture model clustering (using PROC TRAJ (Jones et al., 2001) were analyzed using PROC MIXED (SAS Institute Inc., 2017) cubic repeated measures models were fit using those groupings. Data grouped via rules were analyzed using Repeated Measures and One-Way ANOVAs and Chi Squares. Additional tests included Fisher’s Exact and Binomial Tests.

Results

Sham controls

The Sham control female rats demonstrated a simple paw withdrawal but no HLB to stimulation of the ipsilateral and contralateral hind paws over the 21-day study period, consistent with our prior results for Sham male rats, previously reported (Dean et al., 2017). They were not included in any subsequent analysis.

Initial data exploration

Potential sexual dimorphisms were identified through graphical data exploration, particularly but not only via spaghetti plots (Fig. 1). Possible sexual dimorphisms were a delay in females’ pain development and higher maximal HLB among females. The decision to longitudinally group data was confirmed by observing HLB over time by sex, as the raw data is not generally following a bell-shaped distribution within each time point; combining all males and all females does not holistically represent the collected data (Fig. 1A, B). Secondarily, a possible temporal relationship was sought between HLB and sympathetic (LF/HF) data, indicated in our previously published data (Dean et al., 2017).

Fig. 1

Hyperalgesia-like behavior (HLB) response rate in male (n = 82) and female rats (n = 75) after SNI.

A) Individual HLB response rates: The spaghetti plot of raw data suggests that females begin to develop hyperalgesia behavior later than males, and that some females attain greater maximal hyperalgesia response rates than males. These potential sexual dimorphisms were selected for further analysis. HLB response rates are jittered to avoid overlap. Baseline (pre-SNI) data are shown on day −1.

B) Grouped HLB response rates: Average HLB responses show a delay in hyperalgesia development in female rats. The shaded regions are standard errors of the plotted means (points), not of the line of best fit. Male and female lines had r = 0.95 and 0.99, respectively. Baseline (pre-SNI) averages are not shown because the focus was on hyperalgesia development after injury Combining males and females into single groups is potentially problematic, because the raw data is not bell-shaped at each time point, as seen in A.

Development of HLB after injury

The latent class clustering method identified five clusters for HLB after SNI in combined sex and only male rats, while four clusters were identified for only female rats (n = 70 males, 77 females) (Fig. 2). These data suggest one or more sexual dimorphisms are present in neuropathic pain development. In both sexes, there is a cluster in which HLB response rate was never >20%, and three similarly shaped, increasing clusters, for which the females generally show higher predicted HLB response rates than males. The fifth male cluster displays a sinusoidal shape. The shapes of the combined sex clusters are similar to those of the male-only clusters; these groups are not a combination of the separate-sex clusters, but formed separately. Unless noted otherwise, only results from analysis using the combined sex clusters were reported.

Fig. 2

Clusters of hyperalgesia-like behavior (HLB) response rate after SNI. Using the latent class clustering method on Moving-Average with lag-1 (MA(1)) data, five clusters were identified for males (left panel) and four for females (center panel). Both sexes show three increasing clusters, but females generally show higher predicted HLB response rates. The fifth, male only, group is sinusoidal. These differences suggest sexual dimorphisms are present. Data for males and females were combined for analytic cluster identification and were expected to be comparable but not necessarily identical to those of similarly shaped separate-sex clusters. For combined sex data (right panel), five clusters were identified that were not simply a combination of the separate male and female clusters. As expected values of each cluster are plotted using MA(1) data, the average of days 1 and 3 is plotted on day 2. By definition, values predicted to be <0 were set to 0. Number of animals per cluster are shown in parentheses.

Using the rules-based clustering method, no statistically significant difference in group membership across sex was identified (p = 0.5197) (Table 1). While the discrete mixture model considers both temporal pattern and magnitude of HLB when forming clusters, the rules-method only considers pattern and requires a greater magnitude of change among higher pain scores. The rules-based clustering method was more sensitive to decreases in HLB response rate/ recovery than the latent class method. Thus, the non-significant difference by sex in the rules-based classification does not contradict the sex difference observed via the latent class method. On the contrary, the composition of the different classification schema was not dissimilar (Table 2), although both cluster number (5 latent classes, 9 rules reduced to 7) and definitions differed, particularly with respect to non-development of HLB.

Table 2

Composition of longitudinal clusters by sex.

Combined Sex Latent Class Clustering			Rules-Based Classification
Cluster(s) Description	Malen (%)	Femalen (%)	Category & Definition	Malen (%)	Femalen (%)
Increasing over time (3 groups)	48 (68.6%)	51 (69.9%)	BH: Increasing over time	43 (61.4%)	51 (69.9%)
Rate always < 20	13 (18.6%)	20 (27.4%)	A: Rate always = 0	6 (8.6%)	6 (8.2%)
Sinusoidal pattern	9 (12.9%)	2 (2.7%)	F: Sinusoidal pattern	5 (7.1%)	3 (4.1%)
			CG, D, E, I: Other patterns	16 (22.9%)	13 (17.8%)

The clustering methods were fundamentally different. There were 5 clusters identified using the latent class method, but 9 rules-based groups, which were combined to 7. Lack of development of hyperalgesia-like behavior was defined differently across methods, as shown.

Because the interest was in post-SNI HLB, baseline (pre-SNI) measurements were not included in these trajectory groupings. As mentioned previously, nearly all rats displayed a 0% response rate at baseline.

Development and magnitude of HLB after nerve injury

Mixed models analysis showed that development of HLB was significantly different by sex (p = 0.0319), and by trajectory group (p = 0.0407) (Fig. 3). Importantly, it showed that the pattern of development of HLB differed substantially over time when considering the interaction between sex and HLB trajectory group (p < 0.0001). Within the same group, development of HLB was delayed in females relative to males, and females reached a magnitude of HLB similar to, or higher than, males.

Fig. 3

Mixed model predictions of hyperalgesia-like behavior (HLB) response rate after SNI using combined sex hyperalgesia development clusters. Sex comparisons within the same cluster show the modeled HLB responses rising later for females than males, and some females attaining greater maximal response rates than males. Different HLB development clusters are denoted by different line type (dots, dashes).

Comparisons using repeated measures ANOVA showed similar findings regarding HLB development, specifically difference over time (p < 0.0021) and by rules-based classification group (p < 0.0001). Sex alone did not show a statistically significant difference (p = 0.0997), but the interaction of sex with time did (p = 0.0164). To understand the timing of pain development across sex, a Fishers’ Exact test was employed to compare when each rat’s HLB response rate first exceeded 20%, a limit selected to be outside of the range of random error (Dean et al., 2017, Gemes et al., 2009). Among males that reached a HLB response rate >20%, 49.0% did so on or before the third post-op day, but 81.5% of females did so on or after day 7 (n = 107, p = 0.0006). The split between day 3 or before vs day 7 or later was determined by observing counts by day by sex.

One-way ANOVA comparing the maximum magnitude of HLB response did not find sex differences (p = 0.3142), and the results were similar (p = 0.3092) for sex interacted with the rules-based groups (after excluding group A, “all zeros,” and groups C, E, G, and I for size). Because the ANOVA is a comparison of means, these findings do not contradict those of the mixed model. Regardless, it is noted that only female rats achieved a maximal response rate of 100%.

Longitudinal clustering of LF/HF power ratio

The latent class clustering method identified 3 clusters for sympathetic tone, as indicated by LF/HF power ratios, over time in combined sex and only female rats after SNI, while 2 clusters were identified for only-male rats (n = 39 males, 38 females) (Fig. 4). In both sexes, the lowest-magnitude cluster is U-shaped, but the females’ LF/HF power ratio values are higher than the males’. Both sexes have a LF/HF power ratio cluster with a magnitude around 40; it is the “high” cluster for males and “lowest” cluster for females. The combined sex clusters’ shapes are similar to those of the female-only clusters; these groups are not a combination of the separate-sex clusters, but formed separately. Unless noted otherwise, only results from analysis using the combined sex clusters were reported. No rules-based clustering method was used for this data. Because interest was in the post-SNI sympathetic tone pattern, baseline (pre-SNI) measurements were not included in the trajectory analysis.

Fig. 4

Clusters of LF/HF power ratio after SNI. Using the latent class clustering method, two clusters were identified for males (left panel) and three for females (center panel). Females’ LF/HF power ratio values are higher than the males. These differences suggest sexual dimorphisms are present. Data for males and females were combined for analytic cluster identification and were expected to be comparable, but not necessarily identical, to those of similarly shaped separate-sex clusters. For combined sex data (right panel), three clusters were identified that were not simply a combination of the separate male and female clusters. Number of animals per cluster are shown in parentheses.

Sympathetic tone after nerve injury

Mixed models showed that the changes in LF/HF power ratio were not different by sex (p = 0.1023) or over time (p > 0.5), but trajectory groups (Fig. 5) did show significant differences (p = 0.0036), probably due to the third group present only among females. Within a trajectory, values were sufficiently similar over time to be considered “flat,” probably partly due to high variability and relatively small sample size (n = 77 after excluding those with 2 consecutive missing values, or missing values on days 1 or 21).

Fig. 5

Mixed model predictions of LF/HF power ratio after SNI using combined sex LF/HF power ratio clusters. Sex comparisons show that all three clusters contain females, but only two contain males. For the latter clusters, females have higher LF/HF power ratios than males.

The change in sympathetic tone from baseline using raw data was also explored (Fig. 6). A Fishers’ Exact test was employed to compare when each rat’s LF/HF power ratio first dropped below its own baseline; among males whose ratio dropped below baseline, 81.6% did so on or before the third post-op day, but 42.9% of females did so on or after day 7 (n = 73, p = 0.0160).

Fig. 6

Average LF/HF power ratio among all males (n = 45) and all females (n = 44) over time after SNI. LF/HF power ratio appears to rise after SNI from baseline levels, and then decline, but temporal changes are statistically insignificant. The shaded regions are standard errors of the plotted means (points), not of the line of best fit. Male and female lines had r = 0.54 and 0.97 respectively. Baseline (pre-SNI) averages are not included in the line of best fit formation but are denoted over time via dotted line for visual reference.

Relationship of HLB to autonomic function

Discrete mixture models’ group memberships were compared across HLB and LF/HF power ratio trajectories via Chi Square and showed a borderline statistically significant relationship (p = 0.0972). More moderate LF/HF power ratios were associated with HLB trajectories of lower magnitude (Fig. 7).

Fig. 7

Associated mixed model predictions of LF/HF power ratio and hyperalgesia-like behavior (HLB) response rate using combined sex clusters. The highest-magnitude LF/HF power ratio cluster to lowest are shown left-to-right, respectively. As shown in Fig. 6, no males are identified in the highest-magnitude LF/HF power ratio cluster. These data indicate that more moderate LF/HF power ratios are associated with lower HLB response rates. A temporal relationship may be indicated between changes in predicted LF/HF power ratio and HLB response rate. Predicted LF/HF power ratios were multiplied by 100, to match the scale of predicted HLB response rate for illustrative purposes.

In all cases, however, statistically insignificant trends in LF/HF power ratio trajectory appear to be temporally associated with significant changes in HLB response trajectories (Fig. 7), possibly supporting a biologically rational relationship between sympathetic tone and HLB. These data are included as observations in this exploratory study and merit further investigation.

Comparisons using repeated measures ANOVA showed LF/HF power ratio differences by sex only (p < 0.0001), with females showing a higher magnitude, but not rules-based HLB classification group (p = 0.8857) or time (p = 0.7017). The differences between the two classification methods’ findings are not a contradiction due to the inherent definitional differences in classification scheme, particularly because an ANOVA method is unlikely to be able to differentiate the moderate baseline LF/HF power ratio from those post-SNI due to its natural rise-and-fall of over time (Fig. 6).

Baseline LF/HF power ratio was also explored in relation to the HLB response rate. A Chi Square showed that when LF/HF power ratio dropped below baseline was related to when HLB response rate first exceeded 20% (n = 52; p = 0.0199); LF/HF power ratio dropped when HLB response rate rose in 57.9% of rats, and an additional 21.0% dropped below baseline before HLB responses rose. Further study is merited as these data suggest a temporal relationship between a drop in sympathetic tone and onset of HLB.

Discussion

The present study highlights an understudied aspect of pain behavior, specifically that grouped data illustrating the development of pain measures can be misleading by masking not only different trajectories but also sex-related differences that could have clinical and physiological significance and, if translated to clinical studies, potential diagnostic and therapeutic implications. The grouped data for males and for females demonstrate a robust increase in HLB that is maintained over the 21 days post-SNI, in agreement with the original description of the SNI model (Decosterd and Woolf, 2000). However, our study went further by employing complementary techniques to analyze the individual temporal pattern of hyperalgesia development, as well as comparing behavior in male and female rats. There were several advantages to the analytic tools employed in this study. Classification by temporal trajectories rather than by data magnitude at a specific time-point enables comparisons to focus directly on sexual dimorphisms, even if they were not present over the entire study period.

Both latent class analysis and the rules-based longitudinal classification methods demonstrate that both sexes may follow one of several different patterns in HLB development after a uniform nerve injury. Despite high inter-individual variability, there was concurrence between the findings of the two classification schema, providing evidence for distinguishable sex differences in the development of neuropathic pain after injury. Discrete mixture models recognized five patterns of HLB development in only males and four in only females, suggesting sexual dimorphism. For analysis, patterns identified using both sexes’ data combined were used. The patterns identified with sexes combined were similar to those most populated among the rules-based classification, with the majority of rats, both male and female, demonstrating increased HLB development over the 21 days post injury. However, the varied patterns of response in both males and females likely contribute to the lack of consensus in the literature that report either no sex differences (Dominguez et al., 2009, Severino et al., 2018) or greater sensitivity of neuropathic pain-related behaviors in female rats (Coyle et al., 1995, DeLeo and Rutkowski, 2000, LaCroix-Fralish et al., 2005, Nicotra et al., 2014) after nerve injury.

Heterogeneity of pain outcomes is recognized in clinical situations, such that different or similar neuropathic etiologies may be followed by highly variable pain conditions or the absence of persisting pain (Ephraim et al., 2005), so the variety of temporal patterns we observed is consistent with the complexity expected from clinical observations. A detailed understanding of individual differences in pain development and its pathogenesis could lead to appropriate interventions and improve outcomes. We and others have previously reported similarly high interindividual variability in pain behavior manifestations in various mononeuropathy models in rats, including the absence of HLB in a subpopulation of animals with anatomically confirmed nerve injury (Cui et al., 2000, Djouhri et al., 2006, Gemes et al., 2009, Hogan et al., 2004) An important finding here is that clustered data clearly reveals a group of animals, up to 19% males and 27% females, that do not develop hyperalgesia after SNI. While this has been reported by others in preclinical models where 29–40% rats did not develop neuropathic syndrome after sciatic or spinal nerve injury (Cui et al., 2000, Dean et al., 2017, Gemes et al., 2009, Kupers et al., 1992), such distinctions have rarely been previously described after nerve trauma, likely due to grouping of all data (Decosterd and Woolf, 2000, DeLeo and Rutkowski, 2000, Dominguez et al., 2009, Kupers et al., 1992, LaCroix-Fralish et al., 2005, Roytta et al., 1999). The reasons for variability in, or lack of, development of pain behavior is currently not known but have been attributed to biological, psychosocial, environmental, and genetic factors in humans that are largely absent or controlled in animal studies (Bushnell et al., 2015, Main, 2013, Mogil, 2012, Paller et al., 2009). The variations identified in the present study are unlikely to be due to genetic variation as all the rats were of the same strain, although epigenetic difference may be present, and pharmacologic and anatomic differences in analgesic mechanisms have been noted for identical strains obtained from different vendors (Clark et al., 1992) The observation of individual differences raises the possibility that there are subtypes of temporal patterns of response to injury that are dictated by underlying differences in pathogenic mechanisms. More complex injury models have shown that the operator is a variable that can affect outcomes (Djouhri et al., 2006), and in consideration of human testing variability all the injuries were performed by a single, experienced individual in the present experiments. We have previously observed variations in rat peripheral nerve anatomy at the lumbosacral level (Rigaud et al., 2008) so it is possible that the extent of injury could vary despite the stereotyped surgery. To further limit variability, one individual was responsible for acclimating rats to the testing apparatus and performing the mechanical stimulation. Finally, although we identified a subgroup that lacks sensory abnormality to a noxious mechanical stimulation paradigm, those rats may have responded to a different stimulus e.g. thermal, cold, or non noxious (Hogan et al., 2004, Jaggi et al., 2011, Roytta et al., 1999) and we acknowledge that our findings regarding the specific measure of mechanical HLB may not be transposable to other behavioral tests. Measurement of HLB rather than mechanical thresholds by von Frey filaments may have contributed to the range of patterns as the former is an integrated response regulated supraspinally compared to the latter that is a simple spinal cord reflex. However, as we note above, this test has clinical relevance, and the experience associated with this specific HLB mode of response motivates aversion in rats (Wu et al., 2010), giving it validity that other induced behaviors may lack.

Both mixed models and Fishers exact test identified a delay in the development of HLB in females compared to males. These data are in agreement with findings from Bourquin et al. (Bourquin et al., 2006) that assessed mechanical allodynia-like behavior after SNI in mice and demonstrated higher response rates in males on day 4, representing a dissociation of the response curves that resolved by day 7. Others have also found dimorphism in the time course of allodynic response after chronic constriction injury with females showing a higher threshold for mechanical allodynia compared to males until 17 days post injury (Vacca et al., 2014).

As seen within the separate discrete mixture model classes, the magnitude of response rate reached in females was similar to, or higher, than that of the males, providing further support for sexual dimorphism. However, when the data for all subjects are combined, the sex difference in severity of HLB is inconclusive. This may be a caveat worth considering in the design of human studies.

While sex hormones are often the primary candidates to account for sex differences, this study did not characterize estrous stage effects, in part due to unnecessary stress imposed on the animals from daily smear tests. We could find no previous studies that examined individual development of mechanical hypersensitivity after SNI in male and/or female rats. The present study takes the initial step in detailing sex differences in the development of HLB after SNI. There are reports that estradiol or estrogen receptor β agonists attenuate neuropathic pain behaviors in both female and male rats after nerve injury (Lee et al., 2018, Piu et al., 2008, Vacca et al., 2016, Xu et al., 2019), but they do not assess the contribution of estrogen to sex differences in neuropathic pain. Future studies would be necessary to examine pre-pubescent and post-estrus cycling animals in addition to estrus cycle stage within females in order to examine the role of hormonal state in development of HLB, without speculation.

The variation in pain trajectories, the delay in hyperalgesia development in females, and the potential higher maximum magnitude of HLB suggest that sex differences may contribute to the heterogeneity and that mechanistic differences may shape the development of neuropathic pain in male and female rats. Largely ignored until recently, progress is being made towards an understanding of mechanistic sex differences associated with pain development (Cooper and Craft, 2018, Juni et al., 2010, Mogil et al., 2003, Sorge et al., 2015).

It is well established that sensory and cardiovascular function are closely linked in central pain networks including the periaqueductal gray (Schlereth and Birklein, 2008). In acute stress conditions, sympatho-sensory integration serves a protective role, to maintain homeostasis with an increased sympathetic output coordinated with elevated nociceptive thresholds (Bruehl and Chung, 2004, Nordin and Fagius, 1995, Sheps et al., 1992). While preliminary, the present study suggests that at some point after injury, this balance may be altered, as it is observed that increases in pain behavior occur at the same time as or after a drop in sympathetic tone below baseline. The present findings warrant further investigation to advance our understanding of integrative processes in pain regulation.

The primary limitation of this study was its sample size, although the two methods produced similar findings, supporting differences in the number and shapes of trajectory paths, and the likelihood of sexual dimorphism in pain development. While the sample was large for an animal study it was somewhat small for the modeling methods used. Thus, while we can confidently say that sexual dimorphisms are present, we cannot say that the groups developed here are universally the best representations of the different temporal interactions among sex, and pain-related behaviors.

This study, which minimized variations in surgical procedures and behavioral evaluation while blinding to the extent possible, provides new insight on variability in HLB after nerve injury. Using two complementary analysis techniques, we unmasked patterns of development of HLB after nerve injury in both male and female rats that are not evident when grouped data are analyzed. The discrete mixture model considers both temporal pattern and magnitude of hyperalgesia behavior when forming clusters, whereas the rules-method considers pattern only. The variability in response patterns following nerve injury is similar to that reported in humans, including a variety of temporal patterns and a lack pain behavior. Although the symptoms and signs of pain development in clinical subjects have been studied in detail, the temporal patterns are unexplored yet highly variable, even within pathogenic categories. The present study suggests that there could be mechanistic heterogeneity yet to be identified to explain this variability in manifestations (McCarthy et al., 2012), as has been suggested for convergent/divergent trajectories in pain development between the sexes (Mogil, 2012). Moving forward, the formation of a translational pain research consortium as proposed by Renthal et al. (Renthal et al., 2021) could prove invaluable in providing the samples and tools necessary to understanding mechanisms underlying the variability in pain development.

These data highlight complexity in development of pain-related behavior that should be taken into account when considering treatment strategies, although replication of this study will be necessary to identify if these findings are robust. Consideration of variations in pain development patterns and dimorphic differences may aid the future development of pain therapeutics, including importantly sex-specific treatments.

CRediT authorship contribution statement

Katherine Sherman: Methodology, Formal analysis, Resources, Writing - review & editing, Visualization. Victoria Woyach: Investigation, Writing - review & editing. James C. Eisenach: Conceptualization, Writing - review & editing. Francis A. Hopp: Conceptualization, Writing - review & editing. Freddy Cao: Investigation, Writing - review & editing. Quinn H. Hogan: Conceptualization, Methodology, Writing - review & editing, Supervision, Project administration, Funding acquisition. Caron Dean: Conceptualization, Methodology, Formal analysis, Investigation, Writing - original draft, Writing - review & editing, Visualization, Supervision, Project administration, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.