Fifty female Sprague-Dawley rats weighing 200-250 g were used in this study. Eight rats that underwent no further experimental procedures were randomly selected to provide a normal control (normal group). The remaining 42 rats were assigned to the model establishment. After surgery, rats were screened, and those presenting self-harm behavior were determined to have neuropathic pain syndrome. Twenty-four rats with neuropathic pain were randomly divided into three groups equal in number: model group, EA group, and sham-EA group. All protocols and procedures employed in this study followed animal care guidelines and were approved by the local committee.

BPA established an animal model of neuropathic pain following Wang et al. (2019). Surgical procedures were performed when adequate anesthesia (40 mg/kg sodium pentobarbital intraperitoneal injection) was achieved. Briefly, an approximate 4 cm incision was made along the dorsal midline using the C7 vertebra as a bone mark. To expose the nerve roots from C5 to T1, the lamina was removed from the right side of the vertebra (C4-T1) under an operating microscope (magnification × 10) through a posterior approach. After being dissociated, both dorsal and ventral rootlets from C5 to T1 were avulsed by a retractor. The same experienced operator performed all these procedures.

Before the intervention, a period of three to four weeks was left for the rats to recover from the surgery. A three-day adaption training was then applied to minimize anxiety for the rats. During the EA intervention, rats were held on a customized platform. In the present study, the parameters of EA were developed in accordance with both literature and the clinical experience of acupuncture experts. Additionally, previous studies have shown that EA on “Jiaji” points relieves some pain conditions, including persistent tissue injury (inflammation), nerve injury (neuropathy), cancer and visceral pain (Chen et al., 2018, 2013; Liu et al., 2016; Xu et al., 2010).

For EA intervention, acupuncture was delivered using four disposable, stainless steel needles (diameter 0.3 mm, length 13 mm, Huatuo, Suzhou, PRC) which were inserted into “Jiaji” on the left side (5 mm horizontally away from the 5th cervical vertebra to the 1st thoracic vertebra, depth 5 mm). EA was applied by an electrical stimulator (Huatuo SDZ-II, Suzhou Medical Appliance Factory, PRC) on two non-adjacent needles. The stimulator delivered dense disperse waves of 2/15 Hz at 1 mA (adjusted to the muscle twitch threshold) for 30 min/day. For sham-EA, four needles were inserted into the non-acupoints besides “Jiaji” on the left side. The intervention was performed five times per week for three months.

Mechanical withdrawal threshold (MWT) tests were performed one month and three months post-intervention to assess changes in pain perception. The timeline of this study is shown in Fig. 1. In consideration of the peripheral sensitization mechanism associated with neuropathic pain, the pain threshold of bilateral intact hindpaws was utilized as the indicator of the overall pain threshold aroused by BPA (Rodriguesfilho et al., 2003). Prior study has consistently demonstrated that the MWT of bilateral hindpaws significantly decreases after BPA when compared with pre-modeling (Wang et al., 2019).

Figure 1.

The timeline of the study. A period of three to four weeks was left for the rats to recover from the surgery before the intervention. Resting-state fMRI and mechanical withdrawal threshold (MWT) tests were conducted after one and three-months of intervention.

For MWT tests, rats were individually held in a transparent plastic box on an elevated mesh grid. After a 30-min acclimatization period, a linearly increasing mechanical force (0.6, 1.4, 2, 4, 6, 8, 10, and 15 g) was applied to the plantar surface of the hindpaw by use of Von Frey monofilaments (Stoelting, IL). A response was considered positive if a quick foot withdrawal response, flinching, or licking was observed. The MWT was determined as the lowest force that triggered a positive response to at least 2 out of 4 repeated stimuli (Chaplan et al., 1994).

Resting-state fMRI data were acquired on a 7T system (horizontal-bore animal MRI scanner, Bruker, Germany) at two-time points: at the end of the first and third months of EA intervention. The present study utilized a single transmit and surface receiving coil to receive the signal. Rats were anesthetized with isoflurane (4% initially for induction, 1.5%-2% for maintenance) and given necessary ventilator support during the scanning. Functional MR images were then acquired using an interleaved single-shot echo planar image (EPI) sequence with the same parameters as reported previously (Wu et al., 2018). Scanning parameters included: repetition time 3000 ms, echo time 20 ms, flip angle 90°, field of view 32 × 32 mm² and slice thickness 0.5 mm.

Data were preprocessed by using an SPM8 toolbox (http://www.fil.ion.ucl.ac.uk/spm/) and Matlab 2014a (www.mathworks.com) software platforms. After converting the raw Digital Imaging and Communications in Medicine (DICOM) data to Neuroimaging Informatics Technology Initiative (NIFTI) images, all images were rescaled by a factor of ten such that the image size approximated the scale of the human brain. This allowed the use of processing algorithms originally developed for human data (Moya et al., 2015; Wu et al., 2018). This procedure only changed the dimension information in the header file and required no interpolation of data. Additionally, none-brain tissue was stripped manually to facilitate image processing. Subsequently, acquired functional images were corrected for differences in slice scan time and spatially realigned with a six parameter (rigid body) transformation to remove movement artifacts. A transformation was then estimated from the mean image to a standard template (Schwarz et al., 2006) and applied to fMRI data for individual Montreal Neurological Institute normalization. Spatial smoothing was performed with a full-width half-maximum Gaussian kernel of twice as the voxel size. Subsequently, linear detrending and band-pass filtering (0.01-0.08 Hz) were performed. Finally, nuisance variables (blood oxygen level dependent signal of the global brain, white matter, and cerebral spinal fluid and motion parameters) were further regressed out (Fox et al., 2009).

A functional connectivity matrix was constructed using a graph theoretical network analysis toolbox for imaging connectomics (GRETNA, http://www.nitrc.org/projects/gretna/) (Wang et al., 2015). Schwarz et al. (2006) were used to divide the whole brain into 96 regions, which were then considered as nodes of the brain networks. For each brain region, time series were acquired by averaging the time course across all voxels within each given region. Pearson's correlations between the mean time series of all node pairs were computed as the edges of a brain functional network. The result was a 96 × 96 symmetric correlation matrix for each rat, which was then converted into a binary matrix at a specified sparsity (Suo et al., 2015). Notably, the number of edges generated by this method divided by the maximum possible number of edges is referred to as the sparsity of the network. To obtain networks with the same number of edges, the sparsity range was set for each correlation matrix. This allowed prominent small-world properties in the brain networks to be observed (Watts and Strogatz, 1998). Accordingly, each absolute matrix was restricted at each sparsity level (a threshold range of 0.05 < sparsity < 0.45 with an interval of 0.01) to generate a set of undirected and unweighted binarized networks. Finally, graph theoretical network analysis was performed using GRETNA.

Statistical analysis was performed with Matlab software. Repeated measures analysis of variance was performed to analyze behavioral data. Post-hoc pairwise comparison was performed to identify significant between-group differences at each time point using a Bonferroni-Dunn test. A P-value < 0.05 was considered to indicate a significant difference between samples. The integrated regional network parameters were used for the Kruskal-Wallis test at each time point. If there was a significant overall difference among groups, further group comparison was performed by a nonparametric permutation test at each time point. P < 0.05 after false discovery rate (FDR) correction was considered significant. To determine whether there was a statistical difference in the proportions of more highly connected hub regions, the Friedman test was performed. Wilcoxon signed-rank tests were performed to assess the difference between either two groups at the same time point. To compare the difference between two-time points in the same group, pairwise signed-rank tests were performed. The threshold for significance was set at P < 0.05.

A significant main effect for group comparisons was found for the mechanical withdrawal latency (MWL) of left hindpaws (F (3, 28) = 22.581, P < 0.05) and right hindpaw (F (3, 28) = 13.100, P < 0.05). After one month of EA intervention, MWL of bilateral hindpaws in the model, EA, and sham-EA groups were significantly lower than that of the normal group (P < 0.05). Alternatively, MWL of the left hindpaws in the EA group was significantly higher than that of the sham-EA group (P < 0.01). For MWL of right hindpaws, there was no significant difference among the model, sham-EA, and EA groups. After the three-months of EA intervention, MWLs of bilateral hindpaws in the model and sham-EA groups were both significantly lower than those of the normal group (P < 0.01), while no significant difference was found between the EA and normal groups for MWL in the left hindpaws. MWLs of bilateral hindpaws in the EA group were significantly higher than those in the model and sham-EA groups (P < 0.01) after three-months of EA intervention (Fig. 2).

Figure 2.

Mechanical withdrawal threshold of bilateral hindpaws of four groups after one and three-months of EA intervention. Data are expressed as mean ± standard error of the mean (SEM) for eight animals in each group. △ P < 0.05, as compared to the normal group; ▲ P < 0.01, as compared to the normal group; ◇ P < 0.01, as compared to the model group; * P < 0.01, as compared to the sham-EA group.

Based on the nodal properties, hubs were defined in functional brain networks in each of the four groups of rats at each of the two tested time-points. The hubs identified in each group at each time point are given in Table 1 and mapped onto the rat brain surface for visualization (Fig. 3). Results showed a similar spatial distribution of hubs among the four groups. Hub regions were predominately located in the limbic system (such as entorhinal cortex, nucleus accumbens shell, piriform cortex, hippocampus, and hypothalamus), parietal lobes (motor cortex), midbrain (ventral tegmental area, mesencephalic region, superior colliculus, and substantia nigra), medial temporal lobe (olfactory tubercle, nucleus accumbens, diagonal band and substantia innominata) and subcortical regions (ventral pallidum and thalamus).

Figure 3.

Hub regions of the functional neural networks of four groups of rats after one and three-months of intervention. More significant spheres show regions with a standard deviation higher than the average value of the nodal parameters.

Significant differences in the proportion of more highly connected hub regions were found across all groups (Friedman’s Chi-square 17.07, df 3 28, P = 0.0079). However, further analysis showed that after the one-month of EA intervention, there was no significant difference among the four groups. After three months of EA intervention, the proportion of more highly connected hub regions in the model group was significantly higher than that of the normal group (P = 0.048), while, conversely, the proportion of the EA group was substantially lower than that of the model group (P = 0.049). No statistical difference was found between the EA and normal groups or between the sham-EA and other groups (Fig. 4).

Figure 4.

Differences in the proportion of more highly connected hub regions among four groups after one and three-months of intervention. △ P < 0.05, as compared to the normal group, ◇ P < 0.05, as compared to the model group.

In addition to the abovementioned network metrics, differences of the integrated regional network parameters were determined for each brain area (Table 2).

At the end of the first month of EA intervention, the right mesencephalic region exhibited a significantly higher degree for the model group than for the normal group, while the right mesencephalic region and left periaqueductal grey showed a substantially lower degree for the EA group as compared to the model group. Betweenness centrality of the left amygdala in the EA group was significantly lower than that in the normal group. Concerning nodal efficiency, no significant difference was found among the four groups of rats after one-month of EA intervention.

At the end of the third-month of EA intervention, betweenness centrality of right olfactory nuclei in the model group was significantly lower than that of the normal group, while betweenness centrality of the left anterior dorsal hippocampus and left pons in the model group was significantly higher than that of the normal group. When compared with the model group, the left anterior dorsal hippocampus exhibited a substantially lower degree of hubness in the EA group. Betweenness centrality of right olfactory nuclei and left amygdala in the EA group were significantly lower than that of the normal group. Only for the left anterior dorsal hippocampus was the nodal efficacy in the model group substantially higher than that of the normal group. The right lateral hypothalamus and left amygdala exhibited significantly lower nodal properties in the EA group as compared to the sham-EA group. No significant difference was found among groups concerning degree.