Resting‐state functional connectivity modulates the BOLD activation induced by nucleus accumbens stimulation in the swine brain

Abstract Introduction While the clinical efficacy of deep brain stimulation (DBS) the treatment of motor‐related symptoms is well established, the mechanism of action of the resulting cognitive and behavioral effects has been elusive. Methods By combining functional magnetic resonance imaging (fMRI) and DBS, we investigated the pattern of blood‐oxygenation‐level‐dependent (BOLD) signal changes induced by stimulating the nucleus accumbens in a large animal model. Results We found that diffused BOLD activation across multiple functional networks, including the prefrontal, limbic, and thalamic regions during the stimulation, resulted in a significant change in inter‐regional functional connectivity. More importantly, the magnitude of the modulation was closely related to the strength of the inter‐regional resting‐state functional connectivity. Conclusions Nucleus accumbens stimulation affects the functional activity in networks that underlie cognition and behavior. Our study provides an insight into the nature of the functional connectivity, which mediates activation effect via brain networks.

However, it is unclear which specific neurobiological mechanism mediates broadly distributed BOLD activation patterns. Since BOLD activation can be found across anatomically heterogeneous regions that may not be directly connected to the stimulation site, an anatomical connection would not likely be the exclusive mechanism for delivering the DBS effect, albeit some anatomical connections have been found in DBS-activated regions from diffusion tensor imaging (DTI) and fiber tracing studies in animal models (Britt et al., 2012;Leh, Ptito, Chakravarty, & Strafella, 2007;Lehéricy et al., 2004).
We investigated the pattern and extent of BOLD activation and changes in functional connectivity between activated brain regions after high frequency stimulation in a healthy pig model. In particular, we observed changes in inter-regional functional connectivity (rsFC), focused on brain regions that evoked a significant BOLD response to NAc-DBS, including executive, limbic, thalamic, and sensorimotor networks. Furthermore, we also conducted rsFC mapping (Biswal, Zerrin Yetkin, Haughton, & Hyde, 1995;Fox, Halko, Eldaief, & Pascual-Leone, 2012) for each subject prior to the DBS implantation surgery in order to collect information on the relationship between the strength of rsFC and BOLD activation, as a potential action mechanism of DBS, which has been overlooked in previous DBS-fMRI studies (Gibson et al., 2016;Min et al., 2012;Paek et al., 2015;Settell et al., 2017). The use of a large animal model in our study is presumably more representative of the human brain anatomy (Van Gompel et al., 2011) than small animal models (Albaugh et al., 2016), and therefore, our findings should be assumed to be general in nature, in terms of understanding the therapeutic effects of human DBS. F I G U R E 1 (a) Schematic diagram depicting the experimental procedure and (b) the placement of electrode in the placement of subjects on the pig brain atlas (Saikali et al., 2010). The red dots and colored lines indicate the estimated position of an individual electrode tip and body, respectively (n = 8  Figure 1a outlines the experimental timeline. After initial sedation and intubation, the animal was delivered to a scanner for imaging for an individual anatomical scan (25 min 36 s). The anatomical scan was followed by a resting-state functional MRI (6 min 30 s). After restingstate functional imaging, electrode implantation surgery was followed (2 hr) and the DBS-fMRI experiment was carried out (6 min 30 s).

| DBS electrode implantation
A quadripolar DBS electrode (Model 3389; Medtronic Inc.) was implanted targeting to NAc of the left hemisphere (Knight et al., 2013;Min et al., 2012). The coordinates (arc, collar, and depth) on a stereotactic frame (Leksell; Elekta Co) were determined by COMPASS planning software (COMPASS International Innovations) based on NAc location identified from individual subject's anatomical brain images (Knight et al., 2013) (Figure 1b). Microdrive (Alpha Omega Co.) was used to guide the implantation of a DBS lead. All of these systems can minimize the displacement of an electrode from a targeted position (<1 mm), as shown in our previous phantom tests (Min et al., 2014).

| DBS stimulation parameters and DBS-fMRI acquisition
Following DBS surgery, each animal underwent functional imaging with simultaneous NAc stimulation. Each stimulation block consisted of a 6-s stimulation train followed by a 60-s off period stimulation ( Figure 1a). The block was repeated five times per scan with the initial baseline period. The total time per scan was 6 min and 30 s. The imaging parameters were the same with those in the resting-state scan were used.
The stimulation parameters were as follows: biphasic and bipolar pulse, voltage = 5 V, pulse frequency = 130 Hz, and pulse duration = 100 μs. These parameters were selected based on results from our previous study wherein we found a robust and reproducible BOLD activation in multiple trials across different subjects (Knight et al., 2013;Min et al., 2012;Paek et al., 2015;Settell et al., 2017).

| fMRI data preprocessing
Resting-state and DBS-fMRI data were processed using the AFNI software (Cox, 1996). Animal physiological data were collected to re- Spike removal, slice-timing correction, and within-subject motion correction with six parameters (three translations and three rotations of x-, y-, and z-axis) were applied. Individual subject's anatomical and functional images were coregistered to the pig brain atlas using a 9-parameter linear registration (3-translation, 3-rotation, and 3-scaling) with the cost function of the Hellinger distance (Mallet et al., 2008), and the alignment of coregistration was then visually checked each time. Prior to statistical analysis, spatial smoothing was applied using 3-mm full-width-half-maximum (FWHM) isotropic Gaussian kernel. A temporal filter was not used in this study, because the physiological artifacts were regressed out, and temporal filtering could influence the temporal synchronization between imaging and stimulation onset (Davey, Grayden, Egan, & Johnston, 2013).

| DBS-fMRI BOLD activation map
The general linear model (GLM) analysis was adopted to detect stimulation-induced BOLD activation and measure the amplitude of signal chance. An individual analysis was carried out first, and a group-level (n = 8) statistical analysis was then performed (3dDeconvole and 3dttest, AFNI). In the group-level analysis, significant BOLD activation was detected with the threshold, set as p < .05 (n = 8, t > 3.49, false discovery rate [FDR] corrected). The coordinates, t-scores, and percent signal changes were summarized. The activation map was created with different color encoding based on the percent change in a signal and overlaid on the pig brain atlas in three view planes (axial, coronal, and sagittal). BOLD time courses for each ROI were generated. We first created a sphere mask (radius = 1.7 mm) and applied it on the center location of each ROI. We then extracted all time courses within the mask and averaged them individually, and generated the group-averaged ROI time courses (n = 8). poststimulation-state matrix. Finally, individual subject's CC matrices were normalized by Fisher's Z-transformation and averaged into a group-level CC matrix (n = 8). To detect significant FC during resting state, one-sample t test was applied. To detect the significant change of FC, two-sample t test was used between resting-and poststimulation-state CC matrix.

| The relationship between resting-state functional connectivity and the BOLD response
Pearson's CC (r) and slope of linear regression (s) were estimated between the resting-state functional connectivity (rsFC) of a given ROI to the NAc and the BOLD response of the ROI. In the voxelwise analysis, the same calculation was conducted: the sphere mask (radius = 1.7 mm) was applied to individual voxel in functional networks, and a voxel-wise averaged time course was obtained. Then, statistics (r and s) between rsFC-to-NAc and BOLD response of a given voxel were estimated. Data points were categorized into five functional networks to plot separately.
The correlation and regression analyses were also applied between inter-ROI rsFC and inter-ROI FC change in poststimulation state. The same sort of analysis was applied to analyze the relationship between co-activation of ROIs and inter-ROI FC change in poststimulation state. In those analyses, the absolute CC value was considered.

| Postsurgical evaluation of DBS electrode placement
The placement of the lead and location of electrode tip was examined by reconstructing an electrode-induced image artifact on individual EPI image volumes (Horn & Kühn, 2015). Intensity thresholding was applied to the axial plane of a single EPI image slice. The signal dropout region induced by the image artifact was isolated, and centroids of the region were marked along with the slice direction.
By interpolating those marks, the lead placement could be reconstructed on individual anatomical image volumes. The location of the electrode tip with contacts 0 and 1 was then estimated and mapped on the pig brain atlas (Figure 1b).

| Comparison of EPI signal intensity in stimulation on and off periods
Echo-planar imaging image volumes, obtained during stimulation on and off periods, were averaged across blocks. Voxel-wise signal intensity was extracted within and adjacent area of the signal dropout region. Data points of stimulation on and off period were plotted separately with spatial distance. The statistical analysis (two-sample t test) between the signal intensities of stimulation on and off blocks was conducted to detect the significance of signal difference in conditions. F I G U R E 3 Group-level BOLD activation map of a 130 hz NAc-DBS (n = 8; t > 3.5; p < .05, false discovery rate [FDR] corrected). Activated regions are shown in (a) coronal, (b) axial, and (c) sagittal views on the pig brain atlas (Saikali et al., 2010), and the t-scores for activation are denoted by colors. NAc-DBS induced BOLD activation in multiple cortical and subcortical brain regions. Note that DBS electrodes were implanted in the left hemisphere for all subjects. For abbreviations, see Appendix A

| BOLD activation induced by NAc-DBS
NAc-DBS evoked significant BOLD activation in multiple cortical and subcortical brain regions, but was limited to the ipsilateral hemisphere ( Figure 3). The highest amplitude of evoked BOLD activation was found in the ipsilateral NAc (signal change: 1.27 ± 0.19%, beta coefficient: 1.6, t = 12.11, p < .05, false discovery rate [FDR] corrected), indicating the current spread near the DBS electrode tip likely evoked a direct neuromodulatory effect (McIntyre, Mori, Sherman, Thakor, & Vitek, 2004). However, it must be noted that BOLD activation was also evident beyond the stimulation locus, that is, prefrontal, cingulate, insular, and sensorimotor cortices (p < .05, FDR corrected) (see Table 1 for the summary of activated clusters).
Thus, the results show NAc-DBS induced both local and distal modulatory effects across multiple functional networks including executive, limbic, thalamic, and sensorimotor networks. Interestingly, such effect was primarily lateralized to the left hemisphere, ipsilateral to the DBS electrode (Figure 3b). Only a few regions in the contralateral hemisphere (pIC, CD, Pu, and premotor cortex) showed significant BOLD activation (p < .05, FDR corrected), demonstrating that the DBS effect may be unilateral. Figure 4 shows group-averaged ROI time courses by individual region. The DBS evoked BOLD response was characterized by an initial negative induction (5 ± 3 s) followed by peaks (25 ± 5 s), which is consistent with previously published hemodynamic responses to visual and sensory stimulation (Ogawa et al., 1992(Ogawa et al., , 1993.
Additionally, daCC showed a broad connectivity to many other brain regions including the prefrontal regions, aIC, and NAc as shown in other resting-state FC studies (Cauda et al., 2011;Greicius, Krasnow, Reiss, & Menon, 2003). Our results show that functional coupling could be observed between anatomically separate brain regions under conditions of isoflurane anesthesia.

| Co-activation of ROIs during stimulation
NAc-DBS evoked temporally synchronized hemodynamic responses in ROIs of the ipsilateral prefrontal cortex (r > 0.3, n = 8, t > 2.37, p < .05) (Figure 5b). Importantly, it should be noted that co-activation was found between regions, wherein their anatomical connections are less clear, such as those between the prefrontal cortex and the thalamus. These results support the fact that distal effect of NAc-DBS cannot be fully explained by a direct inter-regional anatomical connection alone, albeit the co-activations in anatomically adjacent regions would not be surprising, that is, regions within limbic system or prefrontal cortex.

| Relationship between resting-state functional connectivity to NAc and amplitude of BOLD activation
We found that a significant positive relationship exists between the rsFC-to-NAc of individual ROIs and their BOLD response (r = 0.52, p < .01) ( Table 1 and Figure 6a). The stronger rsFC to the stimulation locus (NAc) that was presented, the higher was the BOLD response. The relationship between the rsFC and BOLD response, however, was observed, not only in regions of activation above the statistical threshold (p < .05). Rather, the voxel-wise analysis for expanded region of interests further revealed that BOLD response was correlated with the connectivity, that is, voxels in the bilateral

| The change in pair-wise FC connectivity during poststimulation state
Significant decreases in FC (p < .05) were found between many regions in the poststimulation period, that is, daCC, dpCC, pIC, DLPFC, and vaTH (Figures 7 and 8); however, some others were enhanced after stimulation, that is, the pairs of VS-mdPFC, mdPFC-Pu, mdTH-pIC, and mdTH-OFC. Additionally, a few pair-wise FCs became even negative after stimulation: pCC-OFC, dpCC-pIC, AM-pIC, and AM-SAC. The results indicate that NAc-DBS, in general, have an inhibitory effect on major inter-ROI networks after acute stimulation; however, the modulation effect of NAc-DBS could also vary depending on the network being considered.

| The relationship between poststimulation FC change and strength of resting-state FC
While it is not surprising that stimulation could alter the functional connectivity between ROIs, further results showed that the extent

| Assessment of image artifacts induced by electric current
While it is known that the metal components of the DBS electrode can cause magnetic susceptibility artifacts as demonstrated by marked signal reductions in MR images near the implantation site (Settell et al., 2017), it is also possible that the electrical current of DBS could induce an additional image artifact, that is, increasing the size of the signal dropout area. The variation in image signal intensity was assessed near the signal dropout region during stimulation on and off periods; however, no statistically significant difference in signal intensity was found between stimulation on and off blocks ( Figure 10). It therefore appears that the impact of electric current on the images is negligible.

| D ISCUSS I ON
NAc-DBS-induced BOLD activation occurs not only near the stimulation locus (NAc), but also in the distal regions across multiple functional networks that include ipsilateral prefrontal, limbic, thalamic, and sensorimotor areas. While the patterns and extent of BOLD activation were consistent with previous results in pig models and in human studies (Gibson et al., 2016;Knight et al., 2013;Rauch et al., 1994), previous studies have not addressed the mechanism involved in the propagating effect NAc-DBS in the global brain. In our study, several ROIs lacking direct anatomical connections to the NAc showed temporally synchronized BOLD activity during stimulation ("co-activation"), indicating Stim.
that an anatomical connection alone cannot explain the diffuse pattern induced in a global brain. We therefore suggest that functional connectivity between stimulation locus and individual ROI and between ROIs likely play a crucial role in facilitating and propagating the DBS effect in distal brain regions.

| Anatomy and functional connectivity of nucleus accumbens (NAc)
Our resting-state analysis showed that NAc has significant resting-state functional connectivity (rsFC) with the medial portion of the prefrontal cortex (OFC and mdPFC), cingulate cortex, insular cortex, and limbic substructures. These results indicate the presence of significant functional coupling between NAc and an executive network, which is not irrelevant to its anatomical position that links the prefrontal cortex and basal ganglia complex (Brog, Salyapongse, Deutch, & Zahm, 1993;Haber, Kim, Mailly, & Calzavara, 2006;Leh et al., 2007). NAc pathways play a key role in the development of reward-guided behavior by associating reward information with motivational and emotional features of sensory inputs (Groenewegen, Wright, Beijer, & Voorn, 1999;O'Doherty, 2004;Reiss et al., 2005), that is, perceptual learning between visual stimuli and esthetic music reward enhanced the strength of FC between NAc and visual cortex (Salimpoor et al., 2013).

| Local neuromodulatory effects of NAc-DBS
Electrical stimulation directly modulates neuronal activity near the area of electrode implantation. Stimulation altered the spike rate of neurons and their firing pattern near stimulation locus (Anderson & Mullins, 2003;Dostrovsky & Lozano, 2002;Hashimoto, Elder, Okun, Patrick, & Vitek, 2003;McConnell, So, Hilliard, Lopomo, & Grill, 2012). Computation modeling simulated the extent of local effect (i.e., volume of tissue activated) (Butson & McIntyre, 2008;McIntyre & Grill, 1999). Our results also show that NAc-DBS evoked the diffuse activation of BOLD at the outside of the implantation locus ( Figure 3). Thus, the findings suggest that the electric field likely extends from the implantation locus to adjacent subcortical tissues of the NAc, presumably through the local neuronal network, although the extent may depend on the stimulation parameters being applied (e.g., voltage and frequency) (Paek et al., 2015;Settell et al., 2017).

| Global neuromodulatory effects of NAc-DBS
The modulation of electrical stimulation could be limited to a subset of afferent neurons (Canteras, Shammah-Lagnado, Silva, & Ricardo, 1990); therefore, propagation of the effect to distal areas would be probabilistic and sparse in time (Chomiak & Hu, 2007; F I G U R E 4 Time courses for the NAc-DBS-induced BOLD response in the ROIs (all subject averaged, n = 8) in (a) ipsilateral and (b) contralateral hemisphere. Six-second stimulation period is denoted by red bars. The error bars indicate the ±1 standard error of the mean (SEM) of the BOLD signal. The location of each ROI is depicted in the upper right corner of individual figures. For the coordinates, see Table  1. NAc-DBS evoked robust BOLD responses, characterized by an initial negative component followed by a peak (5 ± 3 s and 25 ± 5 s after stimulation onset, respectively). For abbreviations, see Appendix A F I G U R E 5 Pair-wise functional connectivity (FC) is shown in the color encoded matrix (left) and illustrated with lines on a schematic brain diagram (right): (a) resting-state functional connectivity (rsFC) and (b) stimulation-state FC ("co-activation"). Colors indicate the group-averaged Pearson's correlation coefficient (r), ranging between −0.1 and 0.5. Significant resting or stimulationstate functional connectivity (p < .05) is presented as colored lines on the brain illustration. Significant resting-state correlations were found between NAc and multiple cortical and subcortical ROIs. During a stimulation period, multiple ROIs evoked highly significantly correlated BOLD signal activity at p < .05 ("coactivation" F I G U R E 6 The relationship between resting-state functional connectivity (rsFC) to the location of the stimulation (rsFC-to-NAc) and BOLD activation in (a) the ROIs, and all other voxels in (b) ipsilateral (c) and contralateral functional networks. The horizontal and vertical axes denote the rsFC (group-averaged Pearson's correlation coefficient) of each ROI or voxel to NAc and the amplitude of BOLD activation (beta coefficient). The slope (s) of the linear regression is shown as a blue dotted line. Asterisks indicate a significant correlation between rsFC-to-NAc and BOLD activation in a given ROI or voxel of network. A significant linear relationship was found in the ROIs and voxels in ipsilateral functional networks. For the abbreviations, see Appendix A suggest that a certain neurobiological mechanism is operative that allows the local modulatory effect to propagate to the global brain.

| The role of functional connectivity in the propagation of NAc-DBS effect
Some activated areas, on the one hand, appear to be anatomically linked to the NAc. For example, an anatomical connection has F I G U R E 7 Pair-wise functional connectivity (FC) during the poststimulation period presented in the color encoded matrix (left) and illustrated with lines on a schematic brain diagram (right). Colors indicate the group-averaged Pearson's correlation coefficient (r), ranging between −0.1 and 0.5. A significant change in FC was found in subcortical ROI pairs (p < .05), shown as with colored lines on the brain illustration. NAc-DBS suppresses functional connectivity between the ROIs and even alter the direction of connectivity from the positive to the negative in some ROI pairs. See for the abbreviations in Appendix A F I G U R E 8 Significant changes in pair-wise functional connectivity across resting, stimulation, and poststimulation states. The ROI pairs were grouped according to internetwork combinations for visualization (executive-executive, executive-limbic, insular-limbic, limbic-limbic, and limbic-others). The horizontal axis indicates three states, and the vertical axis indicates the functional connectivity (r). The poststimulation correlation coefficient matrix is shown on the upper right corner for reference. Pair-wise FC is increased during the stimulation period ("co-activation") and decreased in poststimulation state ("inhibition"). See for the abbreviations in Appendix A been identified between the NAc and prefrontal regions, as shown by DTI investigations (Britt et al., 2012;Brunenberg et al., 2012;Leh et al., 2007;Lehéricy et al., 2004). Some networks have been identified as involving the basal ganglia complex, that is, an "indirect" (Damoiseaux & Greicius, 2009;Montaron, Deniau, Menetrey, Glowinski, & Thierry, 1996), "hyper-direct" pathway (Brunenberg et al., 2012;Nambu, Tokuno, & Takada, 2002), or the recurrent loops (Leblois, Boraud, Meissner, Bergman, & Hansel, 2006). On the other hand, however, the anatomical connection cannot fully explain the highly spread pattern of activation across large-scale brain networks. Here, our results showed that the stronger the rsFC of a given region to the NAc was, the larger was the induced BOLD response ( Figure 6). Thus, a functional connection that originated from the stimulation locus is highly likely to be the source of the BOLD activation pattern observed in the global brain. Furthermore, a greater FC change was observed in the poststimulation period when the rsFC in the two regions was stronger, supporting the idea that the rsFC may propagate the modulation effect in distal regions from the stimulation locus.

| Functional connectivity change after NAc-DBS
Many inter-regional functional connectivity (FC) changes decrease after stimulation (Figure 7), indicating that the NAc-DBS tends to have an inhibitory effect on functional connections. However, the magnitudes of reduction were not equal in the ROIs, and in some cases, the direction of functional coupling became reversed (i.e., negative correlation) or increased compared with the baseline condition ( Figure 8).
These results suggest that the DBS effect varies depending on the network being considered (Albaugh et al., 2016). It is generally thought that the DBS exerts its therapeutic effect by regularizing or normalizing a pathological oscillation in large-scale networks (Chiken & Nambu, 2014;McIntyre & Hahn, 2010;Rosenbaum et al., 2014;Vitek, 2002).
In line with the "network effect" perspective, we suggest that a variety of effects (inhibition or facilitation of functional connectivity) of DBS Average intensity of stim-on blocks (n=5) Average intensity of stim-off blocks (n=6) 1.7 mm (1 voxel) may reflect its regularizing efficacy for inter-regional neuronal transmission, rather than simply inhibiting network activity.
The current steering technique with a novel lead design (Lehto et al., 2017;Martens et al., 2011) enables the shape of the electric field to be adjusted by the injection of an asymmetric current, thus providing better control of the DBS-induced electric field (Butson & McIntyre, 2008

| The influence of anesthesia on stimulationassociated BOLD activity
The level of sedation (Liu, Zhu, Zhang, & Chen, 2013) and type of anesthetic agent (Angenstein, Kammerer, & Scheich, 2009;Angenstein, Krautwald, & Scheich, 2010) can have an influence on the BOLD hemodynamic response. Higher isoflurane concentrations (the dose over 1.8%) tend to suppress BOLD activity and results in a decrease in spatial specificity with different functional networks being dissolved during resting states (Liu et al., 2013). A burst suppression of the neuronal population was found in cases where high doses of isoflurane were used (Liu, Zhu, Zhang, & Chen, 2010;Vincent et al., 2007), suggesting that isoflurane as an anesthetic agent may have an impact on both hemodynamic and neuronal activity. In contrast to the high concentration, BOLD activation was preserved when a relatively lower dose (<1.3%) was used (Knight et al., 2013;Min et al., 2012;Paek et al., 2015), suggesting that the impact of isoflurane anesthesia would be less influential or negligible when a lower dose is used. It should be noted here that we cannot completely rule out the possibility that BOLD activation and the rsFC measurements were underestimated due to anesthesia. Nevertheless, we observed a robust amplitude (0.3%-1.3% signal change from the baseline) and temporal pattern of BOLD activation that are consistent with previous fMRI results ( Knight et al., 2013;Min et al., 2012;Paek et al., 2015;Settell et al., 2017).

| Clinical implications for human NAc-DBS and Applications for DBS-fMRI
Dysfunctional functional connectivity between the NAc and various brain regions has been implicated in a number of neuropsychiatric disorders (Da Cunha et al., 2015). Combining NAc-DBS with functional imaging allows causal relationships between the stimulation and modulation of functional networks to be identified, thus providing insights into the underlying mechanisms responsible for the therapeutic and adverse effects of NAc-DBS, particulary for neuropsychiatric disorders, for example, OCD and mood disorders (Bewernick et al., 2010;Bewernick, Kayser, Sturm, & Schlaepfer, 2012). Although our subjects were healthy animals, a large animal model can be implicated as a model for the human brain (Van Gompel et al., 2011), providing insights into how the modulation of these networks might underlie the potential therapeutic efficacy of human NAc-DBS.

ACK N OWLED G M ENTS
This work was supported by the National Institutes of Health (NIH C06 RR018898 and NIH R01 NS70872), Hanyang University (HY-2019), and the Grainger Foundation.

CO N FLI C T O F I NTE R E S T
The author(s) declare that there were no conflicts of interest with respect to the authorship or the publication of this article.

AUTH O R CO NTR I B UTI O N S
All of the authors contributed to the study design. Data collection

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Data for: Resting-state functional connectivity modulates the BOLD activation induced by nucleus accumbens stimulation in the swine brain at http://dx.doi.org/10.17632/ r6f82 ypfvt.3 (Cho, 2019).