Altered functional connectivity associated with time discounting in chronic pain

Chronic pain (CP) is a global problem extensively associated with an unhealthy lifestyle. Time discounting (TD), a tendency to assign less value to future gains than to present gains, is an indicator of the unhealthy behaviors. While, recent neuroimaging studies implied overlapping neuro mechanisms underlying CP and TD, little is known about the specific relationship between CP and TD in behavior or neuroscience. As such, we investigated the association of TD with behavioral measures in CP and resting-state brain functional network in both CP patients and healthy subjects. Behaviorally, TD showed a significant correlation with meaningfulness in healthy subjects, whereas TD in patients only correlated with pain intensity. We identified a specific network including medial and dorsolateral prefrontal cortex (PFC) in default mode network (DMN) associated with TD in healthy subjects that showed significant indirect mediation effect of meaningfulness on TD. In contrast, TD in patients was correlated with functional connectivity between dorsolateral PFC (DLPFC) and temporal lobe that mediated the effect of pain intensity on TD in patients. These results imply that TD is modulated by pain intensity in CP patients, and the brain function associated to TD is shifted from a medial to lateral representation within the frontal regions.

. Time discounting in healthy and patients. (a) Plots display individual time discounting functions for healthy subjects (left) and patients (right). Blue and red line represents mean for healthy and patients respectively. (b) Box plot shows log-transformed discounting factor 'k' for healthy (blue) and patients (red). There was no significant difference between healthy and patients (p = 0.74). (c) Scatter plots show the correlation between behavioral measures and time discounting. In healthy subjects (left), time discounting show significant correlation with meaningfulness. In contrast, time discounting only showed correlation with pain intensity in patients (right). R = correlation coefficient.
www.nature.com/scientificreports www.nature.com/scientificreports/ = + SV t kt ( ) 1 1 (1) where t is the delay in days from the reference time point, and k (>0) is the time discount factor. For large values of 'k', subjective value drops rapidly with time; whereas when 'k' is small, subjective value exhibits a slower response. In other words, subjects with small 'k' are more patient and less impulsive than subjects with larger 'k' values ( Supplementary Fig. 1). Behavioral questionnaires of anxiety and depression, insomnia, and fear for movement recorded significantly higher scores in CP patients compared to healthy subjects as previously reported 24,41,42 ( Table 1). Estimated discount factor 'k' (Fig. 1b), 3 subscales of SOC questionnaire (comprehensibility, manageability, meaningfulness), and financial strain had no significant differences between the 2 groups. However, generalized regression analyses revealed that only meaningfulness significantly correlated with TD in healthy subjects, whereas, in patients, only pain intensity significantly correlated with TD ( Fig. 1c, Supplementary Tables 1 and 2).

Brain functional differences between CP patients and healthy.
We assessed functional connectivity on the predetermined 333 cortical parcels derived from boundaries of resting-state functional connectivity (RSFC) 43 and 10 subcortical regions extracted from the Harvard-Oxford Atlas [44][45][46][47] (Fig. 2a). The 343 regions on interest (ROIs) had been assigned to 13 predetermined community networks; auditory (Au), cingulo opercular (CiO), cingulo parietal (CiP), default mode (DM), dorsal attention (DA), front parietal (FP), retrosplenial temporal (RT), salience (Sa), sensory motor (SM) of mouth, sensory motor (SM) of hand, ventral attention (VA), visual (Vi), and subcortical (SC). For the robustness of our brain functional data, standardized graph theoretical method was demonstrated in CP patients, in our healthy subjects, and in 95 age-and gender-matched off-site healthy subjects, taken from NITRC database as reported previously 38 (Fig. 2b). Clustering coefficient, global efficiency, and small-world-ness showed similar trends in the 3 groups, but clustering coefficient significantly decreased and global efficiency increased in CP patients compared to healthy subjects according to each link density, which is a threshold based on percentage of links with higher correlation (Fig. 2c).
Involvement of the default mode network to time discounting in healthy. ROI-based connectivity analysis identified 45 links involving 51 ROIs related to TD, in which each RSFC cluster met a criterion of at least 2 links connected with each other and every link showed less than 0.05 false discovery rate-adjusted p-value (pFDR) (Fig. 3a). When we categorized each link to the predetermined 13 communities, the greatest number of links (6 links of 8 ROIs) were categorized to the DMN. The average of the Fisher's z-transformed correlation coefficients (zr value) of those 6 links are significantly correlated with log (k) (Fig. 3b, Supplementary Table 3). Number of links (degree), clustering coefficient, and efficiency were averaged in the identified 8 ROIs. Repeated measure ANOVA resulted in significant higher degree and efficiency, and lower clustering coefficient in patients compared to healthy (Fig. 3c). However, none of them were associated with log (k) nor meaningfulness across all link densities (Supplementary Table 4). Mediation analysis resulted in a significant indirect effect of meaningfulness on TD through the identified 6 links within the DMN (Fig. 3d, Supplementary Table 5). Averaged zr value was computed within each community and between any two communities. Increased log (k) was significantly associated with increased zr value within the whole DMN including 41 ROIs and decreased zr value between the VA network (23 ROIs) and the Vi network (39 ROIs) (Fig. 3e,f). The differences of graph metrics between healthy and patients are same in the whole DMN as in the TD-related DMN (Fig. 3g). ≥ 2 links), and no overlap was observed with the links in healthy subjects (Fig. 4a). We focused on the right dorsolateral prefrontal cortex (DLPFC) included in the CO network, because the greatest number of connections (7 links) were in the discovered network region. The DLPFC mainly connected to the temporal lobe (TL) negatively, which covered the bilateral Parahippocampal gyrus (Supplementary Table 6). The increased zr value and decreased degree significantly associated with increased log (k) (Fig. 4b, Supplementary Table 7). Moreover, patients showed significant higher degree and lower clustering coefficient compared to healthy (Fig. 4c). Mediation analysis resulted in significant indirect effect of the pain intensity on TD through the coupling of DLPFC to TL (Fig. 4d, Supplementary Table 8). Furthermore, pain intensity was associated with the degree of the whole DMN (Fig. 4e, Supplementary Table 7). Regarding involvement of anxiety and depression in the TD factor and brain functions associated with pain intensity, multiple regression analyses were performed including Kessler Psychological Distress Scale (K6), and model fitting values were compared to the model of pain intensity (pain model). Any models with K6 for log (k), zr value of DLPFC-TL, and degree of the whole DMN showed decreased F-value and adjusted R 2 , increased small-sample-size corrected version of Akaike information criterion (AICc) and Bayesian information criterion (BIC) compared to the pain model (Supplementary Table 9). Finally, we hypothesized that DMN worked on a 3-step mediation model from pain intensity to TD and tested it. The model represents a series of mediation pathway composed of a few brain networks ( Supplementary Fig. 3). The TD-related DMN identified in healthy and the whole connections within DMN were put into the box "A" and www.nature.com/scientificreports www.nature.com/scientificreports/ "B" in the models respectively, as well as every connection within DMN were tested one by one. However, none of them showed significant fitting on these models (Supplementary Table 10).

Discussion
Our findings identified a significant correlation of pain intensity to TD in CP patients. Patients with severe pain showed high discount rate for future value. This finding implies that severe pain will indicate low expectancy for future condition. Eventually, patients with intractable chronic pain have unfavorable cognitions regarding the probability of future experiences, events, and behavior 48 . Therefore, pain intensity is supposed to increase TD and influence an individual's decision making of inter-temporal choices, even if they do not directly relate to pain. In contrast, meaningfulness, a subscale of SOC, correlated with TD in healthy subjects, whereas it was not replicated in patients. This suggests that pain is a more powerful influencer of TD than meaningfulness in chronic pain state. Since CP state shows functional changes of brain neurons, neuro transmitters, receptors, and neural circuits, especially in the brain reward system 49,50 , this neural plasticity may contribute to the different mechanisms of the intertemporal decision-making between healthy and CP patients.
TD variability in patients might contribute to compliance of pain treatments. Higher TD patients will prefer standard pain treatments, such as drug administrations and injection therapy, which can provide their effects within from a few hours to a few days. On the other hand, advanced interdisciplinary pain programs like a cognitive behavioral therapy (CBT) are probably not acceptable for them, because they will need at least a few weeks Only the connection within DMN and between VA and Vi networks showed the significant correlation to the log-transformed TD factor. (f) Scatter plots show association of TD with the whole connection within DMN and the network between VA and Vi. (g) Graph metrics of the whole DMN. Patients showed significantly higher degree and efficiency, and lower correlation coefficient through 2% to 10% link densities. zr = Fisher's z-transformed regression coefficient. All statistical analyses were controlled with age and gender.
www.nature.com/scientificreports www.nature.com/scientificreports/ as a window period to recognize pain modulation. Since negative RSFC of the DLPFC-TL was associated with lower TD in the current study, the negative coupling of the DLPFC-TL may be a potential predictor of successful pain programs. Furthermore, amending higher TD may be an initial therapeutic target for interdisciplinary pain treatment and the DLPFC-TL may be helpful as a biomarker.
DMN including MPFC and PCC is a brain network shown to be active in resting-state 51,52 . In healthy subjects, community-based analysis revealed that TD correlated with RSFC of the whole DMN. MPFC and PCC were brain regions activated by assessment of subjective value 36 . Increased RSFC within a module including MPFC, DLPFC, and PCC was associated with high impulsivity 53 , which is a major personality shaping TD 54,55 . ROI-based analysis identified that TD-associated links in the DMN were composed of DLPFC and MPFC (Supplementary Table 3); that is consistent with previous findings 56,57 . In addition, the increased RSFC of the specific links within DMN represented a negative correlation to meaningfulness. According to a previous study, increased RSFC of the DMN is associated with lower levels of happiness 58 . Furthermore, our result of mediation analysis showed the specific DMN is a mediator on the effect of meaningfulness on TD.
Similar to the behavioral results of the TD-related measures in CP patients compared to healthy, we identified completely different connections associated with TD in patients. Both DLPFC and TL are common brain areas involved in decision making and reward processing 37 . The coupling of DLPFC-TL, however, has never been identified as a TD-related region in healthy. As many previous studies revealed, chronic pain dramatically alters the functional architecture of brain 38,59-66 . Indeed, graph theoretic analyses revealed significant differences of global and local brain function between patients and healthy in this study (Figs 2c, 3c,g, 4c). Furthermore, the mediation effect of the DLPFC-TL on the relationship between pain intensity and TD implies that the functional change of the brain precipitated by CP consequently makes the DLPFC-TL acquire the function related to TD. In other words, as degree of the whole DMN was associated with pain intensity, chronic painful condition makes the impulsive system involved in TD shift from DMN to the coupling of the DLPFC-TL.
Previous literatures identified higher TD in depressive people 67,68 . Since depression is also a confounding factor of pain 69,70 , it may involve to the association of pain to the increased TD and the TD-related brain functions. However, any multiple regression models including K6 did not show improvement of the models' fitting in terms of TD factor and TD-associated brain functions (Supplementary Table 9). This means that anxiety and depression are less effective on TD than pain intensity in CP patients.
Small number of subjects may cause a selection bias for the current study. Unfortunately, this is the first study regarding TD in CP patients so that there is no brain study to discuss the reliability of our findings. However, it can be expected to be a stable result, as our findings in healthy were consisted with previous evidences. www.nature.com/scientificreports www.nature.com/scientificreports/ Since TD data in this study was gathered by a monetary questionnaire, magnitude effect, an amount-dependent discounting, was not controlled. Magnitude effect is a phenomenon that small rewards are discounted more than large ones, and the anomaly is commonly reported in literature 71 . Subjective differences in feeling about the amount of money offered might have influenced outcomes in this study. Financial strain is an objective indicator of socioeconomic status 72 and significantly correlates with income 73 . In the present study, TD factor was indifferent between high and low financial strain people both in healthy (mean, standard error = -1.4, 0.3 vs. -1.5, 0.2; in high vs. low strain) and in CP patients (−1.3, 0.3 vs. −1.7, 0.2). That suggests the magnitude effect has little influence on this study.
Pain intensity likely to make the TD-related brain network shift from the DMN to the DLPFC-TL in CP patients. Although we hypothesized the functional alteration of the DMN precipitated by CP affected the new construction of the TD-related network, the 3-step mediation analyses did not show any significant results ( Supplementary Fig. 3 and Table 7). From the transition of subacute to chronic pain, RSFC between Hippocampus and MPFC dramatically decreases 74 , and Parahippocampal gyrus in TL mediates the relationship between Hippocampus and chronic pain 60 . These brain interactions may be hypothesized to induce the functional shift of the neuro mechanisms underlying TD, and a subsequent longitudinal study will be required to identify that.
In conclusion, CP patients have different resting state brain networks underlying TD compared to healthy subjects. Specifically, in CP patients, the functional coupling of the DLPFC-TL precipitated by chronic pain plays a significant role in TD.

Subjects.
We recruited 19 patients with chronic neck pain (6 males, and 13 females, mean age = 45 years, age range = 21-64 years) and 19 age-and gender-matched healthy volunteers. Specific inclusion criteria for chronic neck pain were (1) experiencing pain in the neck, (2) pain persisted for longer than the last 3 months, and (3) pain intensity of 4/10 or greater on the numerical rating scale (NRS). All participants were right-handed and underwent structural and resting-state functional magnetic resonance image (rs-fMRI) scanning. Subjects were excluded if they reported (1) history of head injury, (2) medication of epilepsy and/or depression, (3) to be unable to keep supine position, (4) contraindication of MRI, and (5) the score of 5/7 or greater on the Stanford Sleepiness Scale (SSS) performed soon after the scanning. All subjects passed radiological diagnosis for brain structural abnormality. In addition, 95 off-site healthy control subjects (30 males and 65 females, mean age = 45 years, age range = 21-65 years) were taken from NITRC 1000 functional connectomes project (http://fcon_1000.projects.nitrc.org/). The present study was approved by the Keio University Institutional Review Board committee (approval number: 20160002) and all participants accepted the written informed consent. All experiments were performed in accordance with Helsinki declaration and an ethical guideline for medical and health research involving human subjects issued by Japanese Ministry of Health, Labour and Welfare. This trial is registered with the University Hospital Medical Information Network (UMIN) clinical trials registry, number UMIN 000024475. Clinical measures. Participants completed self-report questionnaire including NRS for pain intensity, Kessler Psychological Distress Scale (K6), Athens Insomnia Scale (AIS), Tampa Scale for Kinesiophobia (TSK), Sense of Coherence (SOC), financial strain, and time discounting questionnaire. K6 is used for scaling anxiety and depression. TSK includes 11 items to measure fear for movement. SOC represents anti-stressful personality composed of 3 subscales, comprehensibility, manageability, and meaningfulness. Financial strain is rated on a 5-point Likert scale ranging from 1 (very tight) to 5 (enough), and subjects who marked "very tight" of "tight" were categorized into the high-strain group. All questionnaires were administered on the day of brain scanning.
Time discounting factor acquisition. Subjects were asked a set of ten questions in the following format 75 : To me, 'receiving $ X today' is equally as good as 'receiving $ Y in __days, ' where X < Y. Subjects needed to wait longer to get the larger amount of money, $Y, and they filled for the longest acceptable delay (t) that makes the two options the same. The actual amounts of X and Y are one of $5, $10, $15, $20 and $25; the ten combinations in total. Discount function of subjective values, SV(t), was estimated to fit the following Eq. (2) for each subject by the hyperbolic discounting model using non-linear least-square (NLS) method in R 3.4.0 software.
We also tested another assumption using exponential model (3), a time consistent model 71 .
ct where c (>0) is an exponential discount factor. However, it returned unacceptable model fitting (R square < 0) in two subjects. In the other 36 subjects, log-transformation of exponential factor c significantly correlated with hyperbolic factor k, and the R squares in both models significantly correlated each other as well ( Supplementary  Fig. 2). The hyperbolic model and the factor 'k' were therefore employed in the present study.
Brain scanning parameters. Structural  Preprocessing. The preprocessing of each subject's time series of fMRI volumes was performed using the FMRIB Expert Analysis Tool (FEAT, www.fmrib.ox.ac.jk/fsl) and encompassed: Discarding the first four volumes; skull extraction using BET; slice time correction; motion correction; and high-pass and low-pass temporal filtering (0.0075 and 0.1 Hz). These included the six parameters obtained by rigid body correction of head motion, the global BOLD signal averaged over all voxels of the brain, signal from a ventricular region of interest, and signal from a region centered in the white matter were removed from the data through linear regression as non-neuronal fluctuations. All preprocessed fMRI data were registered into standard MNI space, and subsequently registered to the 343 ROIs as mentioned in the main manuscript.
Functional connectivity analysis. Functional correlation maps and connectivity networks were produced using a well-validated method 76 . We first extracted the BOLD time series from a predetermined functional region of interest (ROI) and then computed the correlation coefficient between its time course and the time variability of all other ROIs. Correlation coefficients were converted to a normal distribution using Fischer's z-transform. General linear model (GLM) was computed repeatedly for all connections with a regressor of log-transformed discount factor 'k' under controlling age and gender, and then false discovery rate (FDR) was generated.
Graph theoretical analysis. As we described previously 38,77 , subject dependent threshold was calculated according to the number of edges to compare the extracted graphs. The nodes of each extracted graph comprised with the ROIs and the edges corresponding to the absolute values of correlation coefficient greater than the threshold. We chose 9 values of threshold, from a conservative threshold corresponding to 2% connection density to a lenient threshold corresponding to 10% link density, where link density is the percentage of edges with respect to the maximum number of possible edges [N × (N − 1)/2]. We computed the 'clustering coefficient -a measure of network segregation -and the 'global efficiency' -a measure of the network integration into a community structure of interconnected modules -for each link density using the brain connectivity toolbox (BCT) 78 . They were defined at the nodal level and then the global average was estimated over all nodes. We also computed the 'small-world-ness' based on the tradeoff between clustering and efficiency for each subject. It was computed as the multiplication of the clustering ratio and the efficiency ratio to the random network, which was generated across 100 repetitions in the same number of edges and nodes in each link density 38 . Differences in topological properties between groups were computed using a repeated measure ANCOVA, with age and gender as covariates of no interest.

Network visualization.
ROIs and functional connections were visualized on a surface rendering of a human brain atlas with the BrainNet Viewer (Xia et al., 2013, http://www.nitrc.org/projects/bnv/) 79 . The sizes of the spheres representing their ROI strength scaled by the degrees, and the number of edges. ROIs were color-coded according to the category of the 13 consensus communities. The width of each functional connection was weighted by correlation coefficient of the BOLD signals. Positive and negative correlations were colored orange and green respectively. Mediation analysis. Mediation analysis was performed with PROCESS 3.0 extension in SPSS 24 80 . Indirect effect was computed in bootstrap method permuted 10000 times. The significant fitting for the hypothetic model was defined as the significant non-zeros of the effects in all pathways including the mediators. Age and gender were controlled in all regression analyses.