2.1.

Participants

For this study, we recruited a total of 120 subjects (84 females, $66\pm 8.3$ years) suffering from knee OA. Participants were considered eligible if they had self-reported unilateral or bilateral knee OA pain, according to the American College of Rheumatology criteria.⁴⁴ Nine subjects did not complete the study. The detailed protocol and enrollment procedures are described in our previous article reporting the primary outcomes of our study.¹⁰

Exclusion criteria included a history of brain surgery, brain tumor, seizure, stroke, or intracranial metal implantation; uncontrolled hypertension, heart failure, or myocardial infarction; alcohol/substance abuse; cognitive impairment; and a history of psychiatric illness. The study protocol was approved by the UTHealth Institutional Review Board prior to commencement and is registered with ClinicalTrials.gov (NCT04016272), and participants signed an informed consent form. Each participant was trained in the use of a portable tDCS device ($1\times 1$ tDCS Soterix Inc., New York, United States) to self-administer a 20-min tDCS session using 2 mA for 3 weeks (Monday to Friday) at home or a private room. The tDCS device consisted of two $5\times 7\mathrm{cm}$ electrodes, with the anode placed over the primary motor cortex M1 (contralateral to the affected knee) and the cathode covering the supraorbital area SO (ipsilateral to the painful knee). Participants were randomly assigned to either an active treatment or a sham group in a 1:1 ratio. For the active stimulation, the participant received a 2-mA constant direct current for 20 min according to the anode-cathode placement described before, whereas the sham group only received the electrical current stimulation for 30 s at the beginning and end of the session, with no direct current stimulation for the rest of the time. All sessions were monitored through a secure video conference platform (e.g., WebEx) to ensure that the proper brain stimulation technique was used and to monitor any adverse events.

2.2.

Pain Experimental Protocol and fNIRS Data Collection

Functional connectivity was investigated by recording fNIRS signals during knee-localized thermal stimulation and, on a subset of participants, punctate pain stimulation. Thermal pain stimulation was applied via a commercially available thermal sensory testing machine (Medoc, Inc.) and delivered through a small thermal probe ($1\times 1$ in.) applied to the affected knee. Additionally, mechanical punctate stimulation was applied using a commercially available handheld probe (North Coast Medical, Inc.). All subjects received thermal pain stimulation, whereas only 89 received both thermal and punctate stimulation. The experimental block design consisted of applying pain stimulation for 20 s in six repetitions interleaved with resting periods of 30 s.

Pain-related brain responses were monitored using fNIRS at baseline (i.e., prior to any tDCS session) and after completing the first, second, and third weeks of treatment. The fNIRS device (LightNIRS, Shimadzu, Japan) consisted of eight light sources (semiconductor lasers at 780, 805, and 830 nm) and eight detectors arranged in a grid-like layout covering the PFC and primary motor (M1) and primary somatosensory (S1) contralateral to the affected knee (Fig. 1). This configuration allowed for measuring the hemodynamic activity with 20 optical channels with a source-detector distance of 30 mm approximately at 13.33 Hz.

Fig. 1

Sensitivity map of the prefrontal region and left hemisphere overlaid onto the Colin27 brain model. Sources, detectors, and channels are depicted as red circles, blue circles, and yellow bar links, respectively. The sensitivity profile was computed and displayed with AtlasViewer.

2.3.

Preprocessing of fNIRS Recordings

Prior to fNIRS data analysis, nine subjects were discarded because they withdrew from the study or only completed the baseline (week 0) session. Initial processing of the fNIRS data included the conversion from raw voltage to optical densities and the identification and correction of motion artifacts using the spline interpolation method.⁴⁵ To identify recordings with poor quality data, we applied a quality control analysis on each optical channel via the QT-NIRS toolbox⁴⁶ using the following parameters: ${\mathrm{SCI}}_{\text{threshold}}=0.8$, ${\mathrm{PSP}}_{\text{threshold}}=0.1$, and ${\mathrm{CQS}}_{\text{threshold}}=75\%$. Briefly, the scalp coupling index (SCI) and power spectrum peak (PSP) quantify the strength of the cardiac pulse detectable in a short-timed fNIRS signal (e.g., 3 to 5 s) in an optical channel, and the channel quality score (CQS) represents the time fraction of the entire recording that must have sufficient quality to consider such an optical channel acceptable. Unlike activation analyses in which each channel is associated with an activation value (e.g., a statistical $p$- or $q$-value), functional connectivity analysis is based on the statistical relationship between pairs of channels. Although pruning low-quality channels may make group-level analysis and interpretation challenging because of heterogenous connectivity matrices across subjects (i.e., networks with a different set of nodes), discarding noisy scans in their entirety improves the quality of the datasets without compromising the feasibility of subsequent analytical steps. For this reason, instead of pruning individual channels with low quality (i.e., $\mathrm{CQS}<75\%$), we opted to remove entire scans with fewer than five acceptable channels, as they contributed less than minimally to the functional connectivity analysis. Potential concerns arising from the inclusion of individual noisy channels are alleviated by the use of a robust correlation method⁴⁷ (described later), which in turn is inspired by the robust regression method called autoregressive (AR) modeling using iteratively reweighted least squares (AR-IRLS).⁴⁸ In addition, the quality metrics used in QT-NIRS have been shown to be good indicators of the performance of AR-IRLS.⁴⁹

After data quality evaluation, we applied band-pass filtering (0.01 to 0.2 Hz) to remove undesired components and then converted the data to oxy-hemoglobin (${\mathrm{HbO}}_2$) and deoxy-hemoglobin (HbR). We then removed the remaining systemic noise by regressing out the mean response across channels for each Hb data type and scan (i.e., global system regression).^50–52 This last step was performed due to the lack of short-separation channels in the fNIRS equipment. Finally, we computed the total hemoglobin (HbT) as the addition of the ${\mathrm{HbO}}_2$ and HbR. The entire processing pipeline is illustrated in Fig. 2. Data quality assessment was performed using QT-NIRS,⁵³ and the rest of preprocessing, functional connectivity, and functional activation analyses were conducted using Brain AnalyzIR⁵⁴ and Homer 2⁵⁵ toolboxes and custom scripts in MATLAB v2021b.

Fig. 2

Processing pipeline of fNIRS data.

2.4.

Functional Connectivity Analysis

The brain connectivity analysis typically includes the extraction and representation of the functional connectivity networks. Extraction of the connectivity networks has been done through time-domain methods (correlation or seed-based analysis) or frequency-domain approaches (partial direct coherence or wavelet-based).⁵⁶ The network representation approaches include graph theory methods, connectivity matrices, and atlas-based correlation maps, among others.⁵⁷ However, both extraction and representation methods have different implications in terms of the type of retrieved information and the statistical analysis procedures. Thus, we chose the extraction and representation methods in function of the change in brain integration given the treatments (active versus sham) and the pain stimuli. From the Hb data, we computed the correlation among fNIRS channels using an AR prewhitening approach, which converts the time-series fNIRS signals to meet the statistical assumptions (independent noise terms and zero autocorrelation).⁴⁷ Then we computed the Pearson correlation coefficients for each pair of channels across Hb types and subjects. The correlation was computed from the blocks of pain stimulation, that is, only considering each of the six 20-s trials of the stimulation periods ($\sim 120\mathrm{s}$). In this sense, the correlation value from channels A and B can be understood as the strength of signal association between those channels while the subject is exposed to the pain stimulation. The group-level analysis was performed over the $R$-coefficients using a linear mixed-effects model considering the group-week interaction for each pain stimulation type. In fNIRS connectivity, the selection of which Hb connectivity network to analyze has been debated. In this case, the similarity between the ${\mathrm{HbO}}_2$ and HbR networks may assist in determining whether to use either ${\mathrm{HbO}}_2$ or HbR (when the similarity is high) for the posterior analysis or present the results combining both Hb types. There exist specific metrics to quantify the similarity of connectivity networks; these include the graph edit distance, structural Hamming distance, and differential symmetry index (DSI).⁵⁸ We opted to use the DSI to determine the similarity between the ${\mathrm{HbO}}_2$ and HbR networks because it removes the mathematical contribution of network density that increases the rate of spurious symmetry levels. Once the Hb connectivity network to be used was determined, we estimated the probability density functions (PDFs) from the distribution of the correlation values. The PDFs were computed for each stimulation type and group across weeks to quantify how different or similar the set of functional connections between the sham and active groups were. We employed a kernel density estimation with Gaussian kernels for the density estimation. The similarity between distributions densities was computed by the Hellinger distance.^59,60 For discrete distributions, the Hellinger distance (also known as Jeffreys distance) between two distributions $P$ and $Q$ is defined as $H(P,Q)=\frac{1}{\sqrt{2}}{\Vert \sqrt{P}-\sqrt{Q}\Vert }_2$, where ${\Vert ·\Vert }_2$ denotes the Euclidean norm. We selected the Hellinger distance because it satisfies the axioms of a metric (identity, symmetry, and triangle equality) and generates values in the [0, 1] range, which simplify its interpretation.

To investigate the changes in connectivity at the regional level, we extracted the correlation values at intra- and interregional levels for PFC and primary motor and sensorimotor areas. The regional connectivity was obtained through a subnetwork analysis. First, we obtained the corresponding $t$ values from the correlation values through the expression $t=\frac{r}{\mathrm{S.E.}}$, where $\mathrm{S.E.}=\sqrt{\frac{1-r^2}{n-2}}$ and $n$ is the sample size. Then the statistical significance was computed and corrected using a Benjamini–Hochberg correction.

2.5.

Statistical Analysis at the Subnetwork Level

Statistically, we analyzed functional connectivity both within and between brain regions in an attempt to gain insight into the role of subnetworks involved in functional organization patterns in pain management. Specifically, we described the regional connectivity using mean correlation and mean node degree of subnetworks. The mean correlation represents the strength of connectivity, whereas the mean node degree represents the total wiring cost of a subnetwork.⁶¹ In functional connectivity networks, the division into subnetworks can be done according to different criteria, for instance, anatomical divisions, channels density coverage, and well-known brain circuitry, among others.^62,63 Here we choose a trade-off between the regions of interest relevant to this study and the channel density coverage of our fNIRS setup, thus exploring connectivity in two regions: the PFC and the lateral region that includes the primary motor cortex (M1) and primary sensorimotor (S1) area. We applied the statistical analysis on both measures according to the following scheme. First, the Shapiro–Wilk test was used to examine the normality assumption of data distribution, and then two sample $t$-tests or Wilcoxon Rank sum tests (when normality was not met) were used to compare the active group with the sham group.

3.1.

Data Quality

Our approach to data quality evaluation (i.e., triaging scans with less than five acceptable channels) resulted in the exclusion of 56 out of 799 recordings (or 7% of the total) from further processing. The resulting CQS averaged across the screened dataset is shown in Fig. 3. Expectedly, optical channels over the PFC had a visibly higher quality ($\sim 95\%$) compared with those probing the primary motor and somatosensory cortical areas ($\sim 60\%$), where hair obtrusion contributed to an overall lower signal quality.

Fig. 3

Resulting data quality assessment using QT-NIRS after triaging excessively noisy scans. Optical channels covering (a) the prefrontal region exhibited higher quality (CQS close to 95%) than channels over (b) the sensory-motor brain area (CQS close to 60%).

3.2.

Time Series

In addition to the quality assessment, the removal of the global systemic signals unveils the hemodynamic changes while discarding the effects of physiology from the scalp when subjects are exposed to pain stimulation. Figure 4 shows the block-averaged ${\mathrm{HbO}}_2$, HbR, and HbT group-level responses.

Fig. 4

Global block-averaged response for sham (a) before and (b) after and active (c) before and (d) after applying the global signal regression to remove the global physiological effects of pain punctate stimulation (punctate). ${\mathrm{HbO}}_2$, HbR, and HbT signals are shown in red, blue, and green, respectively.

After removing the global systemic components, hemodynamic changes were expectedly observed in the transition from resting to stimulation periods. Increases in HbT were observed on channels covering the motor and somatosensory regions, whereas negative HbT changes took place on channels covering the PFC. In other words, the signals from channels covering the PFC exhibited an anticorrelated temporal profile in contrast to those from the primary motor and somatosensory regions.

Figure 5 shows the block-averaged hemodynamic temporal profile at the group level for the thermal and punctate pain stimuli at the beginning of the treatment. Different cases of activation and negative activation patterns arose from the pain-evoked responses on both thermal and punctate stimuli. Activation responses are defined as an ${\mathrm{HbO}}_2$ increase concomitant to an HbR decrease pattern, and negative activation responses show the opposite behavior, that is, an ${\mathrm{HbO}}_2$ decrease simultaneous with an HbR increase.

Fig. 5

Block-averaged group-level response. The average hemodynamic response of pain stimulation during baseline (week 1). Thermal stimuli: (a) sham and (b) active; punctate stimuli (c) sham and (d) active. ${\mathrm{HbO}}_2$, HbR, and HbT signals are shown in red, blue, and green, respectively.

With regards to the thermal stimulation type, a considerable set of channels on the primary motor and somatosensory regions showed activation patterns (sham group: S8-D8, S8-D6, S4-D6, S6-D4, and S2-D2; active group: S4-D6, S4-D2, S4-D4, S2-D2, and S6-D4), whereas channels over the PFC exhibited a negative activation response (sham group: S5-D5, S5-D3, S3-D5, and S1-D3; active group: S5-D5, S3-D3, and S1-D3).

With respect to the punctate stimulation, the sham group exhibited more activated channels covering the primary motor and somatosensory regions than the active group (sham: S8-D6, S4-D4, S2-D2, S6-D4, and S2-D4; active: S8-D6, S4-D2, S8-D8, and S2-D2). A similar number of negative activation cases were found on both the sham (S7-D5, S5-D5, and S5-S7) and active (S5-D7, S5-D5, and S5-D3) groups over the PFC region. In contrast, we observed a set of channels with opposite responses to the previous ones, that is, with an activation profile on the PFC and a negative activation response on the primary motor and sensorimotor regions. However, the number of those cases was small, and it was also investigated during the functional connectivity analysis. In general, the activation versus negative activation cases suggested an anticorrelated relationship between PFC and the primary motor and sensorimotor areas that was quantified by connectivity analysis.

3.3.

Functional Connectivity

Functional connectivity among channels was quantified by computing a correlation between the ${\mathrm{HbO}}_2$, HbR, and HbT time series considering the blocks of stimulation periods. At the scan level, we obtained 743 connectivity networks, and at the group level we computed a set of 48 connectivity networks (#stimuli $x$ #groups $x$ #weeks $x$ #Hb types; $2\times 2\times 4\times 3$). The DSI metric was calculated considering only their significant connections ($q\text{-}\text{value}<0.05$). Then we based our analysis only on the HbT results due to the low ${\mathrm{HbO}}_2-\mathrm{HbR}$ symmetry value (mean = 0.02 and std = 0.01). Figure S3 in the Supplementary Material shows the DSI results for all connectivity networks at the scan level, and Figs. S1 and S2 in the Supplementary Material present the group-level HbT connectivity maps (across stimuli, groups, and weeks) and contrasts maps (week 2 to week 0), respectively.

Because we focused our analysis on the variation of functional connectivity strength, we present the functional connectivity models as connectivity matrices. In a connectivity matrix, rows and columns represent the fNIRS channels and the elements $c_{i,j}$ represent the $R$ coefficient value between the channels $i$ and $j$. Note that this representation does not disentangle the directionality of the relationship, but only its strength. Because we did not focus on the directionality of the existing relationships, all of the resulting connectivity matrices are symmetric.

3.4.

Thermal Stimulation

In Fig. 6, we present the HbT connectivity matrices derived from group-level analysis. For each matrix, we show the connectivity matrix and the estimated probability densities. We observed a high-positive intraconnectivity in the PFC and M1_S1 regions in every session for both active and sham groups, as depicted in the first and second rows in Fig. 6 and Fig. S1 in the Supplementary Material. Notably, we found a negative between-regions connectivity. The positive correlation among neighboring channels denotes a localized hemodynamic response, suggesting a modularity organization of the investigated brain regions, whereas the anticorrelated connections exhibit the cooperation between the PFC and primary motor and sensorimotor regions. In addition to the regional organization analysis, we computed the similarity of the connectivity models between groups based on the distributions of the $R$-values computed from each model. Table 1 lists the $H$-distances describing the longitudinal changes across weeks. The minimum distance between the densities was observed at week 0 with $H=0.0250$, whereas the maximum separation was observed in week 2 ($H=0.1254$) (Fig. 6, bottom row). The low $H$-distance found at baseline denotes that the intra- and interregional connectivity in both groups exhibits high similar connectivity strength in response to the thermal stimulation at the beginning of the treatment. On the other hand, the maximum $H$-distance was found in week 2, suggesting the time point with the maximum effect of tDCS treatment. Additionally, important changes can be observed in the density shapes. For the sham group during week 2, the shape shows a more bimodal profile due to a decrease in the density around correlation values near zero, which indicates that some of the weak connections improved their strength. In terms of the network organization, it may be caused by the increasing modularity pattern in the M1_S1 area from week 0 to week 2 [Fig. 6 (second row, M1_S1 area, week 0 versus week 2); Fig. S1 (first row, column 3) in the Supplementary Material]. Additionally, we observed an increase in the interregional PFC-M1_S1 negative connections from week 0 to week 2 [Fig. S1 (first row, columns 1 and 3) in the Supplementary Material]. This increment is confirmed by the statistical contrast map from week 2 to week 0 shown in Fig. S2 (first row, column 1) in the Supplementary Material. On the contrary, the opposite effects were observed for the active group. The increase in the density around zero correlation values suggests that some of the functional links weaken their correlation. This reduced effect in connectivity could be mostly observed in the intraregional connectivity of PFC [Fig. 6 (first row, columns 1 and 3); Fig. S1 (second row, columns 1 and 3) in the Supplementary Material]. The decrease in PFC connectivity was also captured by the contrast map shown in Fig. S2 (first row, column 2, PFC area) in the Supplementary Material.

Fig. 6

Group-level connectivity models of thermal stimulation. Color-coded correlation values for active (first row) and sham (second row) groups across sessions (columns 1 to 4). PFC, M1_S1, and PFC-M1_S1 regions are indicated by the blue, yellow, and green rectangles. The estimated probability densities of the active and sham groups are shown across sessions (bottom row).

Table 1

Hellinger distances between sham and active correlation distributions for the HbT connectivity networks. H(A,B)=0 expresses identical distributions. Maximum differences between distributions across weeks are indicated in bold.

To investigate the interregional connectivity, we applied a threshold of $|R|\ge 0.3$ (where $|·|$ is the absolute value) to the statistically significant functional links ($q\text{-}\text{value}<0.05$). We selected the value of 0.3 because it effectively removes the weak connections in the intraregional connectivity while keeping the interregion links. Figure S4 in the Supplementary Material illustrates the $R$-values of the connections along with the statistical significance before and after applying the $|R|\ge 0.3$ threshold. Then we obtained the mean correlation and node degree within and between the PFC and M1_S1 regions. Table 2 summarizes the statistical comparison between the sham and active groups for every connectivity model at the subnetwork level. At the beginning of the treatment, no statistically significant differences were found. We found significant differences between the sham and active groups in both week 1 and week 2. Only the interregion connectivity network (PFC-M1_S1) exhibited differences between groups during week 1 in terms of correlation. On the other hand, significant differences were found in the PFC intra- and interconnectivity (PFC-M1_S1) subnetworks for both correlation and degree metrics. No significant differences between the active and sham groups in correlation or degree metrics were found in the last week of the treatment.

Table 2

Comparison between groups across weeks for thermal stimulation. Correlation and degree values computed from subnetwork analysis. Statistically significant (p < 0.05) correlation and degree values are indicated in bold.

3.5.

Mechanical Punctate Stimulation

As in the thermal stimulation, a clear positive and negative statistical association was observed in the intra- and interregional brain connectivity, respectively. In general, the M1 and S1 areas exhibited more organized correlated responses for both the sham and active groups than the associations in the PFC (Fig. 7, first and second rows). Additionally, a decrease in the positive correlations was observed in the PFC area for the active group between week 0 and week 2 (Fig. 7, first row, columns 1 and 3), confirmed by their statistical contrast (week 2 to week 0) shown in Fig. S2 (last row columns 1 and 3) and Fig. S2 (second row column 2) in the Supplementary Material. Moreover, an increment in the negative interregional PFC-M1_S1 connections was observed in the sham group between weeks 0 and 2 (Fig. 7 second row, columns 1 and 3) and was also captured by their statistical contrast [Fig. S1 (third row, columns 1 and 3); Fig. S2 (second row, column 1) in the Supplementary Material]. In terms of the density of the correlations, a high overlap is observed between the sham and active groups at the baseline (Fig. 7, last row, first column) with a Hellinger distance of 0.0366. For week 1, the separation between groups increased for the positive and negative correlations; however, in week 2, the distance seemed to be larger for the negative associations but not for the positive correlations (Fig. 7, last row, second and third columns). Finally, the distance between densities started to decrease between groups during week 3. The evolution of the correlations from the sham group showed that, in week 2, the peaks for positive correlations became more evident, suggesting a decrease in weak positive connections (close to zero correlations). Also from week 2, the density of positive correlations for both the sham and active groups started to overlap again, as at the beginning of the treatment. This suggests that, for the punctate pain stimulation, one important difference between groups was captured by the interregional brain connectivity, which exhibited mostly negative correlations.

Fig. 7

Group-level connectivity models of mechanical punctate stimulation. Color-coded correlation values for active (first row) and sham (second row) groups across sessions (columns 1 to 4). PFC, M1_S1, and PFC-M1_S1 regions are indicated by the blue, yellow, and green rectangles. The estimated probability densities of the active and sham groups are shown across sessions (bottom row).

Table 3 presents the comparison between groups across weeks for the punctate pain stimulation at the subnetwork level (correlation and degree values). We applied the same thresholding procedure ($|R|\ge 0.3$) as in the subnetwork analysis. No statistical differences were found at the baseline between the sham and active groups for any subnetwork. However, we found statistically significant differences between the sham and active groups during week 1 for PFC intraregion and PFC-M1_S1 interregion connectivity in the correlation and degree domains.

Table 3

Comparison between groups across weeks for thermal stimulation. Correlation and degree values computed from subnetwork analysis. Statistically significant (p < 0.05) correlation and degree values are indicated in bold.