In our study, 30 LDH patients and 30 healthy controls (HCs) were recruited from the First Clinical Medical College of Zhejiang Traditional Chinese Medical University from 24th November 2020, to 17th August, 2021. Both groups were matched with age, years of education and gender. All the patients were required to complete two clinical assessments, Visual Analogue Scale (VAS) (Hawker et al., 2011), and the Chinese Short Form Oswestry Disability Index Questionnaire (C-SFODI) (Smeets et al., 2011), to evaluate the degree of pain and daily functional activities of participants.

The criteria for inclusion of patients with LDH were as follows: (1) right-handed; (2) aged between 20 to 60 years; (3) evidence of compression of the spinal canal on a lumbar MRI; (4) radiating pain from the lumbar region to the buttocks and lower limb; (5) positive in straight leg-raising test and augmentation test; knee and Achilles jerk reflexes weakened or missing; (6) VAS score ≥ 3/10; (7) C-SFODI score ≥ 20%; (8) not taking pain therapy for at least 1 month before the enrollment.The criteria for inclusion of HCs were as follows: (1) right-handed; (2) aged between 20 to 60 years; (3) no history of LDH; (4) no history of pain caused by any disease and haven’t accepted any pain-related treatment at least 1 month before enrollment.

The criteria for exclusion of all participants were as follows: (1) a history of spinal surgery or severe spinal trauma; (2) bone tuberculosis, tumor, severe osteoporosis, and other orthopaedics diseases; (3) combined with serious medical or psychiatric diseases, such as cardiovascular and cerebrovascular diseases, and those occur in the blood system and digestive system; (4) pregnancy or breastfeeding; (5) having autoimmune diseases, allergic diseases, acute and chronic infectious diseases; (6) having contraindications to functional magnetic resonance imaging (fMRI), such as metal implants, claustrophobia or devices in the body; (7) fMRI examination showing free nucleus pulposus and cauda equina syndrome; (8) vision loss and vestibular dysfunction;(9) acute or chronic pain caused by other diseases.

This study was approved by the Medical Research Ethics Committee and the Institutional Review Board of the First Clinical Medical College of the Zhejiang Traditional Chinese Medical University (No.2017-k-237-01) and was registered in Clinical Trial Registry (No. NCT03475095). We obtained written informed consents from all participants, and the study was conducted in accordance with the principles of the Declaration of Helsinki.

We performed Tuina interventions using traditional Chinese Tuina which includes multiple manipulations (Lee et al., 2017). In the present study, the Tuina interventions were conducted to relax the lumbar muscles by rolling (Figure 1A), kneading (Figure 1B), pushing (Figure 1C) and the pulling and rotating manipulations (Figure 1D) were utilized to correct disordered spinal joints, to relieve pain and to improve lumbar function. Each Tuina session lasted approximately 25 min. All patients received a total of six Tuina sessions, and three times of Tuina proceeded in one week.

Requirements for physician: Tuina was operated by the same physician who has more than 5 years of massage clinical practice experience. Moreover, the physician has received the standardized requirements of this specialty massage training and holds a physician’s certificate. Before conducting treatment for patients with LDH, the physician had to standardize the specific Tuina operation according to the operation essentials to maintain stability in strength, shape and frequency.

Requirements for Tuina environment and equipment: Multi-functional massage beds were used, along with proper adjustment of angle and height which were set to keep patients comfortable. Meanwhile, standard treatment towels and disposable bed sheets provided by Zhejiang Provincial Hospital of Chinese Medicine were used. In our study, all patients enrolled completed six Tuina sessions, and had no complications from Tuina sessions.

All participants underwent a 3T Siemens scanner (Verio, Siemens AG, Erlangen, Germany) with a 12-channel head coil. HCs group underwent one fMRI scan, while the LDH group needed to be scanned twice, before Tuina and after the sixth Tuina. The functional images were obtained using a gradient echo-planar imaging sequence. The scanning parameters were as follows: repetition time (TR) = 2,000 ms, echo time (TE) = 30 ms, 43 interleaved axial slices, matrix size = 64 × 64, field of view (FOV) = 220 mm × 220 mm, flip angle (FA) = 90°, slice thickness = 3.2 mm, gap = 0 (voxel size 3.4 × 3.4 × 3.2), number of volumes = 230. The structural images were obtained by a sagittal T1-weighted 3D sequence with magnetization prepared rapid gradient echo (MPRAGE). The scanning parameters were as follows: TR/TE = 8100/3.1 ms, sagittal slices, slice number = 176, matrix size = 256 × 256, FOV = 256 mm × 256mm, FA = 8°, slice thickness = 1 mm, gap = 0 (isotropic voxel size = 1 × 1 × 1). During rs-fMRI scanning, all participants were asked to keep their eyes closed, not to think about anything and not to fall asleep.

The data preprocessing were performed using RESTplus V1.25 software (Jia et al., 2019) and included the following steps: (1) removing the first 10 volumes; (2) slice-time correction; (3) realignment; (4) normalization was performed using T1 image new segment; (5) removing the linear trend; (6) nuisance regression, including the white matter, the cerebrospinal fluid and Friston-24 head motion parameters (Friston et al., 1996); (7) filtering in the frequency range between 0.01 and 0.08 Hz; (8) spatial smoothing with a 6-mm full width at half maximum Gaussian smooth kernel (only for static and dynamic functional connectivity analysis).

The Kendall’s coefficient of concordance (KCC) was used to measure the local synchronization of the time series of neighboring voxels as follows (Zang et al., 2004):

where W is the KCC among given voxels, ranging from 0 to 1; R_i is the sum rank of the ith time point; R¯ (n+1) K)/2 is the mean of the R_i; K is the number of time series within a measured cluster (K = 7, 19, and 27, respectively. K = 27 is used in the current study); n is the number of ranks. The ReHo value of each voxel was then divided by the global mean ReHo of each participant for standardization purposes. Noted that the spatial smoothing with a 6-mm isotropic full width at half maximum (FWHM) Gaussian kernel was performed after ReHo calculation.

The significant clusters detected by ReHo analysis between LDH-pre patients and HCs were selected as seed regions of interest (ROIs) to calculate static functional connectivity. Briefly, for each participant, the mean time course of each ROI was correlated with the time courses of each voxel using Pearson correlation. The Fisher’s z-transform was applied to convert FC maps into Z maps to acquire optimum normality.

The significant clusters detected by ReHo analysis between LDH-pre patients and HCs were selected as seed regions of interest (ROIs) to calculate dynamic functional connectivity (dFC). The sliding window approach was applied to obtain the dFC maps for each participant by DynamicBC (V2.2)¹ (Liao et al., 2014). We selected a window length of 50 TRs (100 s) and a window overlap of 98% (step size by 1 TRs) to compute the dFC of each participant (Li et al., 2018, Ma et al., 2021). The FC map for each participant were computed within each window, generating a series of FC maps. Subsequently, the coefficient of variance (CV) of dFC maps across time was calculated to measure the temporal variability of intrinsic brain activity. Finally, the dFC variability of all participants were then transformed into standardized z scores by subtracting the mean and dividing by the SD across each voxel to enhance data normality.

The demographics and clinical variables were analyzed by the Statistical Package for the Social Sciences (SPSS) 22.0. The differences between the LDH patients and the HCs in age, years of education and clinical scores were tested with Student’s t-test. The gender difference was tested using the Pearson Chi-Square test.

For the ReHo analysis, we performed the two-sample t-tests to compare LDH-pre and HCs. Gaussian Random Field (GRF) was applied in multiple comparison corrections (voxel-level p < 0.05, cluster-level p < 0.05). The significant whole clusters obtained from ReHo analysis were defined as ROIs. For the static FC and dFC analysis, two-sample t-tests was also applied to compare LDH-pre patients and HCs. GRF was applied in multiple comparison corrections (voxel-level p < 0.05, cluster-level p < 0.05).

To evaluate the Tuina effect, the abnormal brain regions identified by group comparisons (LDH-pre vs. HCs) were created as brain masks separately. The mean ReHo and FC values (both static and dynamic) were extracted within the brain masks. Two-sample t-tests were performed to compare these values between LDH-pre and HCs, and those between LDH patients after Tuina interventions (LDH-pos) and HCs (p < 0.05). Paired t-test was used to compare mean ReHo values between LDH-pos and LDH-pre (p < 0.05). GraphPad Prism 8 was used to assess the changes of ReHo/FC values.

In total, 3 LDH patients and 2 HCs who had maximum head movement exceeding 3 mm or 3°were excluded from the subsequent statistical analysis. we ultimately included 27 LDH patients (17 males, age: 32.2 ± 9.5) and 28 HCs (17 males, age: 31.8 ± 8.1). The two groups had no significant differences in age (p = 0.8552), years of education (p = 0.5691) and gender (p = 0.8638). Compared with the LDH-pre patients, the LDH-pos patients had significantly lower VAS scores (p < 0.0001) and C-SFODI scores (p < 0.0001) (Details are shown in Table 1).

The result of ReHo analysis showed that LDH-pre patients had decreased ReHo values in the left orbital part middle frontal gyrus (LO-MFG, Montreal Neurological Institute (MNI) coordinate: −30 48 −6) (Table 2 and Figure 2).

LO-MFG was selected as a seed to compare the FC between LDH-pre and HCs. We didn’t find any significant difference regions between LDH-pre patients and HCs after GRF multiple comparisons correction.

LO-MFG was selected as a seed to compare the dFC variance between LDH-pre and HCs. Two-sample t-test revealed LDH-pre patients shown decreased dFC variance in left Fusiform, increased dFC variance located in left orbital inferior frontal gyrus and left precuneus, compared to the HCs (Table 3 and Figure 3).

We extracted mReHo values within LO-MFG mask from LDH-pre, LDH-pos and HCs separately. LDH-pos patients had increased ReHo compared with LDH-pre patients (p = 0.029908). Meanwhile, both the LDH-pre and LDH-pos patients showed significantly decreased ReHo compared with the HCs (p < 0.0001 and p < 0.001) (Table 4 and Figure 2).

We extracted mean dFC variance values within left Fusiform, left orbital inferior frontal gyrus and left precuneus mask from LDH-pre, LDH-pos and HCs separately. In left Fusiform, LDH-pos patients had increased dFC variance compared with LDH-pre patients (p = 0.009), no significant difference compared with HCs (p = 0.9796). Meanwhile, LDH-pre showed significantly decreased dFC variance compared with the HCs (p = 0.0004). In left orbital inferior frontal gyrus, LDH-pre patients had increased dFC variance compared with HCs (p < 0.0001), no significant difference compared with LDH-pos and HCs (p = 0.0827 and p = 0.1891). In left precuneus, LDH-pre patients had increased dFC variance compared with HCs (p = 0.0002), no significant difference compared with LDH-pos and HCs (p = 0.1279 and p = 0.0529) (Table 5 and Figure 3).

In this study, rs-fMRI was used in combination with ReHo and seed-based FC analysis to examine the coherence of local neural activity and whole brain connectivity in LDH patients. Next, using this information as a basis, we evaluated the effective of Tuina intervention. Compared with HCs, patients with LDH had reduced ReHo in LO-MFG and anomalous variance of dFC between LO-MFG and three brain regions, including the central brain region of the default mode network (the left precuneus), the medial prefrontal cortex (mPFC) (the left orbital inferior frontal gyrus), and the left fusiform cortex. After Tuina treatment, the coherence of LO-MFG neural activity of LDH patients increased and their dFC variability was also improved.

We finded that he ReHo in the LO-MFG, located in the mPFC, was lower in LDH patients, agreeing with the result from Zhou et al. (Zhou et al., 2019). The reduced ReHo in LO-MFG could indicate that this gene is involved in the pathology of LDH and that its ability to regulate emotion is deteriorating (Zhao et al., 2019). The decreased activity of the LO-MFG may result in a weakened impact on mood regulation, an increase in anxiety and depression, and a reduced ability to use emotion adjustment techniques to alter the amygdala’s response to harmful stimuli (Rolls, 2019).

The mPFC form a hub (or core) node of the DMN, and disrupted connectivity of the DMN in chronic low back pain (Huang et al., 2019; Ng et al., 2021). We selected the LO-MFG seed to calculate static and dynamic functional connectivity. The dFC variance of left orbital inferior frontal gyrus and left precuneus were larger in LDH patients compared to HC. The left fusiform dFC variance was lower in LDH patients than in HCs. The fusiform, a portion of the temporal and occipital lobes, is known to be engaged in several neuronal pathways that are involved in recognition (Jonas et al., 2015; Beelen et al., 2022), memory (Liu et al., 2021; Yan et al., 2022), and emotion (Liu et al., 2021; Yan et al., 2022), despite the fact that its role is not entirely understood. The orbital inferior frontal gyrus and the precuneus are both a member of the standard mode network. When the brain is at rest, the standard network serves as the main hub of functional connections across brain regions, monitoring both internal and external environments, episodic memory, and ongoing cognitive and emotional activity (Brewer et al., 2011). Numerous studies have shown that individuals with chronic low back pain have aberrant brain areas and functional connections associated to miswiring (Zhang S. et al., 2014; Ng et al., 2021). Additionally, Zhang and colleagues demonstrated that patients with cLBP had increased ALFF in the hippocampal/parahippocampal gyrus but a decrease in ALFF values in the remaining gyrus DMN regions when their spontaneous low back pain increased after the pain amelioration maneuver (Zhang et al., 2019). In comparison to healthy controls, individuals with discogenic leg pain had decreased ALFF values in DMN areas, according to Zhou and colleagues (Zhou et al., 2018). This means that a DMN problem may also be the root of the cognitive and behavioral issues identified in LDH.

Tuina has a strong analgesic response when used to treat chronic LDH (Kong et al., 2012). Previous research on the mechanism of analgesia has mainly concentrated on the biomechanical effect, which is hypothesized to be the main therapeutic mechanism. For example, Du et al. investigated the effect of lumbar manipulations with finite element technique, showing the relative displacement between the intervertebral disc and adjacent nerve roots and the internal stress change of the lumbar intervertebral disc (Du et al., 2016). But for pain relief, the central nervous system must be mediated (Ng et al., 2021). The default mode network is part of the triple network model which explains the whole range of pain related co-morbidities (De Ridder et al., 2022). To clear the effect of Tuina on the default mode network in LDH patients is of great importance to understand the central mechanism of Tuina for pain relief.

In the past, the resting magnetic resonance studies on the effect of Tuina mainly focused on changes in spontaneous brain activity or functional connectivity, but not on whether Tuina can restore abnormal brain activity or functional connectivity to a normal level. Actually, to prove that Tuina can reverse the brain function damaged by disease, is of great significance to reveal the remodeling of nervous functions mechanism of Tuina effect. We finded that after effective Tuina intervention, the ReHo value of the LO-MFG decrease in LDH patients increased. The dFC variability of the default mode network anomalies in LDH patients also improved at the same time. Our results lend credence to the notion that the DMN activity is the neural correlate of LDH and that the control of dFC in important DMN brain areas may help to explain the analgesic effect of Tuina and the brain response of LDH patients to pain management after Tuina.

Our study has several drawbacks as well. First, because we applied tight inclusion and exclusion criteria, our sample size was small. Future studies should examine local synchronous activity anomalies in LDH patients using a larger sample size. Second, since there was no intervention in HCs, we assumed that the fMRI data has no change within two weeks of the study period, but this may not be the case. In future research, a group of HCs receiving manipulation of intervention can be added. Last, following the sixth TMS, we did not obtain any follow-up information from patients with LDH. Future studies might involve lengthier follow-ups to obtain data on the duration of pain relief after Tuina.