Transcranial Direct Current Stimulation Modulates Connectivity of Left Dorsolateral Prefrontal Cortex with Distributed Cortical Networks

Transcranial direct current stimulation (tDCS) is widely used as a neuromodulation tool for clinical and research purposes (for reviews, see Dedoncker, Baeken, De Raedt, & Vanderhasselt, 2021; Filmer, Mattingley, & Dux, 2020; Razza et al., 2020; Galli, Vadillo, Sirota, Feurra, & Medvedeva, 2019; Osório & Brunoni, 2019; Bennabi & Haffen, 2018; Lucchiari, Sala, & Vanutelli, 2018; Strobach & Antonenko, 2017; Dedoncker, Brunoni, Baeken, & Vanderhasselt, 2016a, 2016b; Jantz, Katz, & Reuter-Lorenz, 2016; Kim, Ekstrom, & Tandon, 2016; Nitsche & Paulus, 2011; Paulus, 2011; Nitsche et al., 2008; Wagner, Valero-Cabre, & Pascual-Leone, 2007; Gandiga, Hummel, & Cohen, 2006). For example, several tDCS studies have targeted left dorsolateral prefrontal cortex (lDLPFC) to modulate long-term memory (Mizrak et al., 2018; Leshikar et al., 2017; Sandrini et al., 2014, 2016; Javadi, Cheng, & Walsh, 2012; Javadi & Walsh, 2012; Zwissler et al., 2014), working memory (Talsma, Broekhuizen, Huisman, & Slagter, 2018; Talsma, Kroese, & Slagter, 2017; Hill, Fitzgerald, & Hoy, 2016; Mancuso, Ilieva, Hamilton, & Farah, 2016; Trumbo et al., 2016; Carvalho et al., 2015; Andrews, Hoy, Enticott, Daskalakis, & Fitzgerald, 2011; Keeser et al., 2011; Fregni et al., 2005; Marshall, Mölle, Siebner, & Born, 2005), cognitive training effects (Au et al., 2016; Martin, Liu, Alonzo, Green, & Loo, 2014), and attention and cognitive control (London & Slagter, 2021; McIntire, McKinley, Nelson, & Goodyear, 2017; Nejati, Salehinejad, Nitsche, Najian, & Javadi, 2020; Nelson, McKinley, Golob, Warm, & Parasuraman, 2014; McKinley et al., 2013; Gladwin, den Uyl, Fregni, & Wiers, 2012). Prefrontal cortex is thought to contribute to these functions by using information about behavioral goals to modulate activity in posterior cortical areas (Miller & Cohen, 2001; Tomita, Ohbayashi, Nakahara, Hasegawa, & Miyashita, 1999). Accordingly, to the extent that stimulating lDLPFC produces behavioral effects, we would expect that stimulation should modulate potential interactions between lDLPFC and other cortical areas. However, our understanding is still sparse on how the tDCS might influence the interregional dynamics of lDLPFC.

We investigated whether tDCS applied to lDLPFC can modulate network-level functional connectivity between lDLPFC and other regions of the brain. We used the temporal correlation of BOLD signal fluctuations as the measure of functional connectivity between brain regions (i.e., resting-state connectivity; Cole, Ito, Bassett, & Schultz, 2016; Cole, Bassett, Power, Braver, & Petersen, 2014; Vincent, Kahn, Snyder, Raichle, & Buckner, 2008; Fox & Raichle, 2007; Damoiseaux et al., 2006; De Luca, Beckmann, De Stefano, Matthews, & Smith, 2006; Vincent et al., 2006). The spatial structure of these correlations suggests that the brain is organized in multiple semimodular “small-world” networks characterized by high within-network connections and sparse between-network connections (Ashourvan, Telesford, Verstynen, Vettel, & Bassett, 2019; Power et al., 2011; Bullmore & Sporns, 2009; Kaiser, Martin, Andras, & Young, 2007; Sporns & Zwi, 2004), and we aim to investigate how stimulation modulates the target region's connectivity within the network structure. Prior studies have shown that tDCS applied to lDLPFC influences general connectivity strength of the target or neighboring networks. In Keeser et al. (2011), participants went through two experimental sessions, one with active stimulation and the other with sham stimulation, and each stimulation session was preceded and followed by fMRI scanning. Authors showed that the stimulation induced significant increases in lDLPFC activation and significant changes of connectivity within networks that included brain regions that are close to lDLPFC. However, the connectivity changes specific to lDLPFC or between lDLPFC and the networks were not directly tested in this study. In another study, connectivity analyses using lDLPFC as the seed region showed that the stimulation increased connectivity between lDLPFC and bilateral parietal regions inferior parietal lobule [IPL], superior parietal lobule [SPL]); however, this study focused on pairwise functional connectivity relationships and did not take the network structure into consideration (Mondino et al., 2020).

We investigated how lDLPFC stimulation modulates the way in which the target region interacts with other parts of the brain network structure. Given that lDLPFC has been identified as a node in the “frontoparietal control network (FPCN)” (Schaefer et al., 2018), we considered the possibility that tDCS might increase the relative importance of lDLPFC within the FPCN compared to its importance across other networks. An alternative possibility is that stimulation might increase the breadth of lDLPFC connectivity across different networks beyond FPCN. To test these predictions, we conducted a combined stimulation and imaging experiment, in which different groups of participants underwent an active or sham stimulation. Functional connectivity between regions was estimated using resting-state fMRI scans acquired prior to and following the stimulation, and used to compute graph theory metrics that characterize a node's connectivity in the brain network structure. Finally, we examined how stimulation altered the network relationships of lDLPFC.

A total of 60 healthy participants completed this study (Day 1 study session; see Procedure section for details). Participants were active duty Air Force military members recruited from Wright-Patterson Air Force Base. Inclusion criteria for this study included the absence of neurological or psychological disorders, head injury (concussion history, traumatic brain injury), recent trauma or hospitalization, impairment of vision, hearing or motor control, dependency on alcohol, caffeine, or nicotine, and any current medication that may affect cognitive functions. Of the 60 participants who completed the session, nine participants were excluded from analysis because of excessive motion during fMRI scanning (> 3 mm) and a total number of 51 participants were included in the final analyses (13 female, 38 male). Participants were randomly assigned to the sham, 1 mA, or 2 mA active stimulation group (N_sham = 22, N_active = 29 [N_1mA = 13, N_2mA = 16]). All data were collected at Dayton Children's Hospital, and all experimental procedures were approved by the Air Force Research Laboratory Institutional Review Board at Wright-Patterson Air Force Base and were fully described to the participants before they consented to participate in the study.

This study was conducted as part of a larger multiday testing study. This study focused on noncumulative, immediate effects of stimulation using data from Day 1 resting-state fMRI. Figure 1 depicts the procedure of the pretesting and Day 1 sessions. One or 2 days prior to the actual testing, participants came in for the consenting procedure and task practice. They were given a 5-min-long practice training on the Mackworth Clock test (McKinley, 2018; McIntire et al., 2017) that was to be used during stimulation. On each testing day, participants went through a stimulation session and two MRI sessions: one preceding and one following the stimulation session. An imaging session included one block of resting-state fMRI, three blocks of task fMRI (n-back task), structural MRI, diffusion tensor imaging, magnetic resonance spectroscopy (MRS), and arterial spin labeling (ASL). During resting-state fMRI scanning, participants were instructed to keep their eyes open and to focus on a fixation dot at the center of the screen. Participants were taken out from the scanner after the first (prestimulation) imaging session and guided into a behavioral testing room for a stimulation session. A scanning session was about 70 min long, and the average interval between the start of the prestimulation and poststimulation scanning sessions was 158 min. Poststimulation resting-state fMRI was conducted within 1 hr from the stimulation, and this interval is well within the duration that our stimulation regime shows behavioral effects for (McIntire, McKinley, Goodyear, & Nelson, 2014). Nevertheless, one may argue that tDCS effects might wear out before the scanning and mere epiphenomenal effects of stimulation (i.e., thinking about the tingling sensation during the stimulation session) could result in false-positive findings. In order to address such possibility, we obtained participants' rating on the arousal level (fatigued vs. energized) and four types of sensation that are often associated with stimulation (see under Stimulation section below for details).

Transcranial DC stimulation was delivered with a MagStim DC stimulator (magstim.com, Whiteland). Instead of the standard wet sponge electrodes, custom, composite electrodes were used as anode and cathode (Sherwood, Madaris, Mullenger, & McKinley, 2018; McIntire et al., 2014, 2017; McKinley et al., 2013; Park, Hong, Kim, Suh, & Im, 2011). As in the prior studies that used composite electrodes, a composite electrode consisted of five Na/NaCl EEG electrodes that were arranged in a circular pattern with an inner diameter of 1.6 cm. EEG electrodes in the array were separated by 0.1 cm from the neighboring ones (measured from outer edge to outer edge). The anode was placed on the lDLPFC (approximately F3 in International 10–20 EEG system), and the cathode was placed on the contralateral (right) bicep. Electrical field modeling using the Finite Element Method has shown that this approach delivers the current evenly distributed among the five EEG electrodes and, when used with a 2 mA current intensity, delivers stimulation with an estimated current density of 0.199 mA/cm² to the target area (McIntire et al., 2017; McKinley et al., 2013; see Figure 2 in McKinley et al., 2013, for the Finite Element Method model of estimated current distribution). Conductive gel was applied on the electrodes to ensure current conduction, and electrodes were secured to the target area using medical bandages. A constant current of 1 mA or 2 mA was delivered for 30 min in the active stimulation group and for 30 sec in the sham stimulation group. Electrodes remained in place for the full 30 min in both groups.

One minute into the stimulation, participants completed a sensation and fatigue questionnaire. On this questionnaire, participants rated the intensity of four types of sensation (heat, itching, pain, and discomfort) on a 10-point scale, where 0 corresponded to feeling nothing and 10 meant unbearable. The protocol was to turn off the stimulation in the case one reports a sensation of 7 or higher on any category, but there was no incident. Fatigue was rated on a 7-point scale, with 1 feeling most fatigued and 7 feeling energized.

For the remaining duration of the 30-min stimulation, participants completed a vigilance task (an adapted version of Mackworth Clock Test; McIntire et al., 2017). This was done as neuronal modulatory effects via electrical stimulation are largely influenced by the neuronal state before and during the stimulation (Silvanto, Muggleton, & Walsh, 2008). Specifically for tDCS, prior research indicated that anodal stimulation applied while the target region is activated by task engagement has increased effectiveness (Filmer, Lyons, Mattingley, & Dux, 2017; Ruf, Fallgatter, & Plewnia, 2017; Andrews et al., 2011). Left DLPFC, our stimulation target region, is highly engaged in attentional/cognitive control functions such as vigilance (Kim, Kim, & Im, 2017; Nelson et al., 2014; Langner & Eickhoff, 2013) and working memory (Jansma, Ramsey, de Zwart, van Gelderen, & Duyn, 2007; Ranganath, Johnson, & D'Esposito, 2003; D'Esposito, Postle, Ballard, & Lease, 1999; D'Esposito et al., 1998; Barch et al., 1997; Braver et al., 1997).

MRI images were collected using a 3T GE Discovery scanner with a 24-channel head/neck coil at Dayton Children's Medical Center. T1-weighted structural images were acquired using a spoiled gradient recalled acquisition sequence (flip angle = 16°, FoV = 25.6 cm, image dimension: 256 × 256 × 164, voxel size = 1 × 1 × 1 mm). Resting-state functional images were acquired using gradient EPI sequences (TR = 2000 msec, TE = 20 msec, flip angle = 90°, FoV = 24 cm, image dimension = 64 × 64 × 41, slice thickness = 3 mm, voxel size = 3.75 × 3.75 × 3.5 mm, top–down interleaved). The functional image acquisition included four initial dummy scans to ensure signal stabilization, and the scanning lasted 12 min and 14 sec. Dummy scans were discarded, and a total number of 363 volumes were further processed for analyses.

Preprocessing and data analyses were performed using SPM12 (www.fil.ion.ucl.ac.uk/spm) and the CONN toolbox (Whitfield-Gabrieli & Nieto-Castanon, 2012). Functional images were slice-time corrected for differences in slice acquisition times, realigned, normalized, spatially smoothed (FWHM = 8 mm), and segmented using a local-global parcellation atlas (Schaefer et al., 2018; the atlas provides parcellations of 10 different levels of granularity ranging from 100 to 1000 parcellations, and we used the atlas with 200 parcels).

Data were denoised using the component-based noise correction (CompCor; Behzadi, Restom, Liau, & Liu, 2007) method and scrubbing as implemented in the CONN toolbox. Specifically, the following time series were used as temporal covariates in the first-level analysis and removed from the BOLD functional data using linear regression: 12 time series of the estimated motion (three translation and three rotation parameters, and their first-order derivatives), BOLD time series from subject-specific white matter (top 5 PCA parameters) and CSF (top 5 PCA parameters), and a binary vector that marked outlier time points identified by ART (www.nitrc.org/projects/artifact_detect). ART outliers were defined as images with a framewise displacement greater than 0.5 mm from the previous image and images with a global mean intensity more than 2 SDs away from the mean image intensity of all scans. These thresholds are slightly more conservative than the default thresholds in the CONN toolbox (0.9 mm, 5 standard deviations), which are considered to be “intermediate” thresholds (97th percentile). On average, approximately 17% of scans were outliers (ranges: 2%–53%). The resulting BOLD time series were band-pass filtered (0.009 Hz < F < 0.08 Hz) and averaged within each ROI. Fisher-transformed correlations between the resulting filtered time series of ROI pairs were used as the final estimates of connectivity between ROI pairs.

Graph theory measures that characterize a node's connectivity strength (node degree) and connectivity distribution (participation coefficient) were computed using the resting-state connectivity measures and an independent, predefined network structure. For valid testing of the stimulation effects on graph measures, it was crucial to obtain graph measures based on a common, objective network structure across individuals and scanning sessions. To achieve this, the following steps were taken: 1) We acquired stable solutions of modular/network structures from our resting-state connectivity data using a consensus clustering approach (see below for details), and 2) we computed graph measures by assessing the outcome connections against a network structure of our ROIs that was predefined using resting-state functional connectivity from a large sample (n = 1489; Schaefer et al., 2018). The Schaeffer atlas assigns ROIs into networks of two different granularity—seven and 17 networks—and we used the seven-network solution that includes visual, somato-motor, dorsal attention, ventral attention, limbic, frontoparietal control, and default networks.

The functional connectivity network structure from each individual and session was obtained using a consensus clustering approach (Lancichinetti & Fortunato, 2012), which allows for obtaining stable community solutions when the community detection method is not deterministic. A consensus matrix was obtained as following. First, an adjacency matrix of ROI-by-ROI connectivity strength (Fisher's Z) was obtained for each participant and session (see Resting-state fMRI Data Processing section). The adjacency matrices were entered into a clustering algorithm (community_louvain) implemented in BrainConnectivityToolbox (sites.google.com/site/bctnent; Rubinov & Sporns, 2011), and this procedure was repeated 1000 times. For each ROI pair, the ratio that the two ROIs were partitioned into the same community out of the 1000 solutions (i.e., “consensus” rate) was computed, and the consensus rates were thresholded at 0.035 (consensus rates smaller than the threshold were set to zero). The outcome ROI-by-ROI matrix of thresholded consensus rate was used as an input to the next iteration. This procedure (clustering an ROI-by-ROI matrix 1000 times, obtaining a consensus matrix, and thresholding the consensus matrix) was repeated iteratively until the resulting matrix reached a complete consensus (all ROI pairs would have a consensus rate of 0 or 1). We considered an “edge” (i.e., connection) to exist between ROIs if the consensus rate was 1 between the pair of ROIs.

Using the predefined network memberships of ROIs and the consensus matrix solutions that are described above, we estimated graph measures that quantify a node's connectivity strength with given networks (node degree) and how well-distributed across different networks a node's connectivity is (participation coefficient). Connectivity strength of an ROI (i.e., a network node), node degree, was estimated as total number of the node's edges divided by the number of all possible edges. For example, Brodmann's area (BA) 9's node degree within FPCN was computed as the number of edges between BA 9 and other nodes within the FPCN divided by the number of possible edges between BA 9 and all FPCN nodes (i.e., the denominator is N-1 in a network with N nodes). Participation coefficient was computed using the participation_coef_sign function in BrainConnectivityToolbox.

Statistical analyses were conducted using R statistical toolbox 3.4.3. (www.r-project.org/) and custom MATLAB (www.mathworks.com) codes. One-way between-group ANOVA was used to test whether participant groups (sham vs. active stimulation) differed in their age, scanning time intervals, or the sensation and fatigue ratings.

A general linear mixed-effects model (GLMM) and permutation-based nonparametric tests were used to test for stimulation effects on node degree and participation coefficient of targets. Graph theory measures (i.e., node degree) were entered as the predicted variable, and stimulation group (sham/active), scanning session (pre-/poststimulation), and a Group × Session interaction were entered as predictor variables, with the subject variable included in the model as a random intercept. To determine whether the variability of graph measures is explained by the stimulation, the statistical significance of the interaction term (Sham/Active × Pre-/Poststimulation interaction) was critically tested using permutation tests. A permutation distribution of interaction coefficients was acquired by running the regression model 1000 times after shuffling group membership (active/sham) and session (pre-/poststimulation) labels. Specifically, each participant was randomly reassigned to one of the groups (sham/active) while keeping the group size the same (N_sham = 22, N_active = 29). Session labels (pre-/poststimulation) were flipped in randomly selected participants. Shuffling variable labels in this manner guaranteed that the within-subject variability was kept intact. The interaction coefficient obtained from unshuffled, real data was tested against the permutation distribution. The sign of the interaction term is irrelevant to the interpretation of results here; therefore, coefficients smaller than the 5th or greater than the 95th percentile of the distribution were considered statistically significant.

Stimulation groups were comparable in age and the intervals between the scanning sessions prior to and following the stimulation (Table 1; all p > .1). On average, both groups reported ratings below 2 (from 0 being nothing and 10 being unbearable) on all four types of sensation, and importantly, the rating did not differ between groups for any type (all p > .1). In addition to sensation, participants also reported how fatigued or energized they felt on a 7-point scale (with 1 being most fatigued and 7 being energized). The fatigue level that the sham and active stimulation groups reported were identical (sham group = 4.09, active group = 4.10; p > .1). Therefore, it can be assured that effects of the stimulation were not driven by participants' awareness of the stimulation condition. When the active stimulation group was further broken into 1- and 2-mA stimulation intensity groups, the three groups (sham, active-1 mA, active-2 mA) did not differ in any metrics reported here (Table 2).

Among the 200 parcels defined in the Schaefer atlas, ones that were centered in BA 9 or BA 46 (MacDonald, Cohen, Stenger, & Carter, 2000) were identified as lDLPFC. This approach identified two neighboring parcels in the left prefrontal cortex—a posterior dorsal one corresponding to BA 9 (violet in Figure 2; Figure 2B for volumetric views) and an anterior ventral one mapping to BA 46 (blue in Figure 2; Figure 2C for volumetric views). Both of these ROIs have been identified as nodes in the FPCN (Schaefer et al., 2018; Figure 2A). In addition, to test the specificity of stimulation effects, we used the rDLPFC as a control region. Three parcels were identified as rDLPFC using the same approach as for the lDLPFC described above, and all three nodes belonged to the FPCN.

We first examined how the stimulation influenced the within-FPCN connectivity strength of lDLPFC ROIs. To test this, GLMM was conducted on lDLPFC nodes' within-FPCN connectivity strength with the stimulation group (active/sham), scanning session (pre-/poststimulation), and an interaction term as predictor variables (see Statistical Analysis section for details of the regression model). GLMM permutation tests revealed that within-FPCN connectivity strength was significantly increased by stimulation for BA 9 (|z| = 1.90; Figure 3A), but not for BA 46 (|z| = 1.39; Figure 3B) or any of the rDLPFC control ROIs (all |z| < 0.7). We additionally investigated whether the BA 9-FPCN connectivity change was stimulation intensity-dependent. Specifically, we tested the effects of the stimulation intensity (Intensity × Session interaction) on BA 9-FPCN connectivity using GLMM permutation tests. The analysis revealed no significant difference in the stimulation-driven BA 9-FPCN connectivity changes between the 1- and 2-mA stimulation groups (interaction term |z| = 0.11).

Next, we assessed stimulation effects on the stimulation target's connectivity characteristics outside FPCN. We approached the outside-network connectivity via two different measures: the connectivity strength and the distribution of connections outside FPCN.

First, outside-FPCN connectivity strength of an lDLPFC ROI was quantified as node degree that only took nodes that are not part of FPCN into consideration. For example, BA 9's outside-FPCN connectivity strength was computed as the number of edges between BA 9 and non-FPCN nodes divided by possible number of edges between BA 9 and all non-FPCN nodes (i.e., denominator was the total number of non-FPCN nodes). GLMM permutation tests revealed that the connectivity strength of the lDLPFC ROIs (BA 9, BA 46) with nodes outside FPCN was not significantly altered by the stimulation (BA 9: |z| = 1.23, Figure 4A; BA 46 |z| = 0.29, Figure 4B). Similarly, stimulation did not alter rDLPFC nodes' connectivity with nodes outside FPCN (all |z| < 0.8). In summary, the stimulation did not alter the net strength of the connectivity between lDLPFC ROIs and brain regions outside FPCN.

In our next analyses, we assessed the effects of tDCS on the participation coefficient, the measure that characterizes how well-distributed a node's connections are (Meunier, Lambiotte, & Bullmore, 2010; Bullmore & Sporns, 2009; Guimerà, Sales-Pardo, & Amaral, 2007; Guimerà & Amaral, 2005) across the seven networks that were predefined in the atlas (Schaefer et al., 2018). GLMM regression revealed that the participation coefficient of BA 46 was reduced by lDLPFC stimulation relative to sham (|z| = 1.98; Figure 5B), whereas stimulation did not significantly affect participation coefficient values in BA 9 (|z| = 0.37; Figure 5A) or the rDLPFC (|z| < 0.2). We also tested whether the modulation of BA 46 participation coefficient was stimulation intensity-dependent, and found no difference in the stimulation effects between the two intensity groups (|z| = 0.25).

Close examination of the data indicated potential between-group differences in the baseline (prestimulation) BA 46 participation coefficient values (Figure 5B, light blue), and we further inspected this result. Consistent with this impression, we found that prestimulation BA 46 participation coefficient values were significantly higher for participants in the active stimulation condition relative to participants in the sham conditions (|z| = 1.82). To address whether the interaction effect was merely a consequence of prestimulation group difference, we conducted follow-up analyses on subgroups of sham and active condition participants that were closely matched in the prestimulation participation coefficient values. Participation coefficient-matched groups were identified as follows: First, we identified the range of the BA 46 participation coefficients that overlapped between the active and sham stimulation groups, and excluded participants that fell outside the range. This eliminated a few participants from each stimulation group, and yielded 20 and 27 participants for the sham and active stimulation groups, respectively. Next, for the participation coefficient value of each participant in the sham stimulation condition, the closest participation coefficient value was identified from active stimulation participants. The search for a matching pair was conducted through the sham participants' participation coefficient values in a random order. This matching process was repeated 1000 times with a new random order each time. The goodness of matching was assessed as the mean difference of the participation coefficient, and the matching solution with the smallest mean participation coefficient difference was selected for further analysis (difference mean = 0.0258; standard deviation = 0.0291). Permutation GLMM revealed that 1) BA 46 participation coefficients were matched between the sham and the active stimulation groups (|z| = 1.26) in the selected subsample and that 2) stimulation significantly reduced the BA 46 participation coefficients (|z| = 1.68, significant interaction effects) in this subset of participants.

The goal of this study was to delineate how tDCS applied to the lDLPFC modulates interactions between the lDLPFC and networks of the brain. We used functional connectivity measures computed from resting-state fMRI scans acquired before and after the participants underwent a stimulation session. Relative to sham stimulation, active stimulation altered the functional connectivity of lDLPFC in two ways. First, stimulation increased the functional connectivity between BA 9 and other nodes within the FPCN. Second, stimulation impacted the connectivity of BA 46 across different networks. While the net connectivity strength of BA 46 was not altered within FPCN or outside FPCN, the connectivity distribution (participation coefficient) of BA 46 was decreased by stimulation. In contrast to the effects of stimulation on lDLPFC, stimulation had no significant effects on network-level functional connectivity metrics in homologous rDLPFC regions that are also part of the FPCN.

Our results build on prior findings showing that lDLPFC stimulation impacts the connectivity of FPCN and the functional connectivity of lDLPFC. Keeser et al. (2011) examined stimulation-induced changes in the functional connectivity of resting-state networks. Using independent component analysis on resting-state fMRI scans obtained from 13 healthy participants, they identified four networks: default mode, left frontoparietal, right frontoparietal, and self-referential networks. Next, authors used dual regression method (Nickerson, Smith, Öngür, & Beckmann, 2017; Beckmann, DeLuca, Devlin, & Smith, 2005) to identify subject-specific spatial maps for the networks, which were entered into the group-level statistical tests to assess stimulation effects in a network (contrast: (after real tDCS > baseline1) > (after sham tDCS > baseline)). Analysis of the left frontoparietal network revealed increased coactivation between regions within the superior and inferior frontal gyri, the inferior parietal lobule, and the posterior cingulate gyrus. Apart from the connectivity changes, authors also found significant increase of local activation in the anodal target region (BA 6). However, the independent component analysis did not identify the stimulation target as part of the left frontoparietal network, and as a consequence, the dual regression analysis did not assess the connectivity relationship between the stimulation target and other regions in the network. In this study, we capitalized on a network atlas that is predefined from a large, independent sample (Schaefer et al., 2018), and tested how stimulation impacts the stimulation target's network relationships with FPCN and other networks.

Whereas Keeser et al. (2011) focused on the connectivity strength changes in different brain networks, other research focused on pairwise connectivity between lDLPFC and other brain regions. Using a seed-based connectivity analysis, Mondino et al. (2020) showed that the lDLPFC stimulation increased resting-state functional connectivity between lDLPFC and bilateral parietal regions (IPL, SPL). IPL and SPL belong to FPCN, and therefore these findings are consistent with our finding that lDLPFC connectivity within FPCN was increased by stimulation. In this study, we further took the network structure into consideration, and found that the lDLPFC stimulation might increase the role of lDLPFC in the FPCN via modulating the connectivity between lDLPFC and other regions within the network.

Given that the participation coefficient reflects the distribution of connectivity across networks, it is somewhat counterintuitive that stimulation significantly strengthened connectivity of BA 9 within FPCN but stimulation did not significantly affect its participation coefficient. We suspected that this might reflect a concomitant increase of connectivity between BA 9 and other networks in addition to the FPCN. The numerical increase of outside-FPCN connectivity of BA 9 (Figure 4A) was consistent with this idea, and when we further probed BA 9's connectivity with each network, tDCS numerically increased connectivity between BA 9 and all other networks with the exception of Salience Ventral Attention Network. In other words, BA 9 obtained stronger affinity/functional connectivity with a broad range of networks, which did not significantly alter the relative distribution of its connections across networks.

In contrast, there was a significant Group × Session interaction effect on the participation coefficient of BA 46. This finding is interesting as it suggests that the stimulation may not only modulate connectivity within the targeted network but also modulate the target's connectivity with other networks. Given that BA 46 connectivity differed between the sham and active stimulation groups in the prestimulation data (see interaction effects and prestimulation differences in Results section), we performed a follow-up analysis on subsamples of participants who were equated in the prestimulation participation coefficient values. Even after matching prestimulation participation coefficient values, we found a significant effect of stimulation on the BA 46 participation coefficient, suggesting that our finding was not simply because of sampling error. The finding that the connectivity distribution (participation coefficient) of BA 46 was decreased by stimulation suggests that BA 46 became less likely to be interacting with multiple other networks.

There are some limitations of this study. The study design included only sham and active stimulation groups and did not have an active control condition. Although the current design with a sham stimulation condition allows for comparing effects of the stimulation to the absence of stimulation, we cannot conclude from this study that lDLPFC is the only stimulation target that can lead to the connectivity changes shown in our results. Relatedly, one may argue that connectivity changes might be mere epiphenomena of stimulation-related sensation, or arousal. We ruled out this possibility via demonstrating that there was no between-group difference in the sensation and arousal ratings. However, we acknowledge that an optimal approach is to include an active control condition in the design and that studies following up on the present findings should include an active control condition. Another caveat is that the active stimulation group included participants who received different stimulation intensities (1 mA or 2 mA). We used these intensities as lDLFPC stimulation delivered at 1 mA or 2 mA has been shown to be effective for yielding behavioral effects in prior research (e.g., for vigilance/attention, 1 mA: Nejati et al., 2020; Nelson et al., 2014; Gladwin et al., 2012; 2 mA: McIntire et al., 2017; McKinley et al., 2013), but we treated them together as the active stimulation group in the analyses as we did not have strong a priori predictions on the effects of stimulation intensity. Nevertheless, to address the possibility that stimulation effects shown in this study could have been driven by stimulation of one intensity and not the other, we first assured that there was no preexisting difference between 1 and 2 mA groups (Table 2; demographic feature, scanning intervals, and stimulation-related sensation reports). More importantly, we tested whether the two groups differed in any of the stimulation effects using GLMM and permutation tests, and found no between-group difference (see Results section for details). Therefore, it is unlikely that the connectivity changes observed in this study were predominantly driven by one stimulation intensity group.

With the present data, we are not able to directly relate the observed network changes to cognitive performance measures and, therefore, cannot draw conclusions on what would be the functional benefit of the lDLPFC-FPCN connectivity changes. Building on the present results, future studies can investigate whether lDLPFC stimulation improves cognition by modulating connectivity of BA 9 within the FPCN. Accumulated fMRI findings have established that the lDLPFC plays essential roles in various cognitive functions including cognitive control (meta-analysis and reviews: Vanderhasselt, De Raedt, & Baeken, 2009; Neumann, von Cramon, & Lohmann, 2008; Nee, Wager, & Jonides, 2007; Owen, McMillan, Laird, & Bullmore, 2005), long-term memory processing (meta-analysis and reviews: Blumenfeld & Ranganath, 2007; Fletcher & Henson, 2001; Nolde, Johnson, & Raye, 1998) and working memory maintenance (meta-analysis and reviews: Wager & Smith, 2003). Therefore, on the one hand, increasing neuronal excitability in the lDLPFC via anodal tDCS might be sufficient for modulating these task functions. On the other hand, it is notable that these cognitive functions also engage other regions in the FPCN (Lemire-Rodger et al., 2019; Roberts, Libby, Inhoff, & Ranganath, 2018; Thakral, Wang, & Rugg, 2017; Spaniol et al., 2009; Blumenfeld & Ranganath, 2006; Owen et al., 2005; Wager & Smith, 2003; Dobbins, Foley, Schacter, & Wagner, 2002; Pessoa, Gutierrez, Bandettini, & Ungerleider, 2002). This brings up the possibility that the tightened lDLPFC-FPCN connectivity may contribute to behavioral modulations associated with lDLPFC stimulations. Future studies using a wider array of cognitive tasks should allow for directly investigating this question.

The present results raise the question as to whether the kinds of modulations seen following lDLPFC stimulation would generalize to targeted stimulation of other nodes of FPCN. For example, if tDCS is applied to a node of FPCN other than lDLPFC, would the connectivity of that node be modulated the same way as lDLPFC? Or is it a unique characteristic of lDLPFC that its network relationships are particularly malleable? Future research should also address the importance of the duration and additivity of modulation effects. Are the modulation effects on the target connectivity long-lasting? If so, what does the decay function look like? In addition to the duration of stimulation effects, it would be informative to understand whether the observed stimulation effects can be cumulative over repeated stimulation. A prior transcranial magnetic stimulation study showed that stimulation repeated over five consecutive days (i.e., five stimulation sessions) brought about memory enhancement that lasted for a long term (e.g., 15 days; Wang & Voss, 2015). Some behavioral and imaging studies have tested the effects of repeated tDCS, but so far evidence is sparse for additive effects. tDCS repeatedly applied with a 10-hr-long gap between stimulation sessions did not have an additive behavioral benefit (McIntire, McKinley, Goodyear, McIntire, & Nelson, 2020). tDCS applied repeatedly over three consecutive days significantly influenced brain perfusion measured using arterial spin labeling (Sherwood et al., 2018); however, we do not know whether these tDCS effects additively increased across sessions in a dose-dependent manner. Future research is needed to understand whether the target connectivity modulation effects can be sustained and/or enhanced via repeated application of tDCS.

We thank Casserly R. Mullenger and Aaron T. Madaris for their assistance in data collection.

Reprint requests should be sent to Kamin Kim or Charan Ranganath, Center for Neuroscience, University of California Davis, 1544 Newton Cir., Davis, CA, or via e-mail: [email protected]; [email protected].

Kamin Kim: Conceptualization; Formal analysis; Investigation; Software; Visualization; Writing—Original draft; Writing—Review & editing. Matthew S. Sherwood: Conceptualization; Data curation; Investigation; Methodology; Project administration; Writing—Review & editing. Lindsey K. McIntire: Data curation; Investigation; Project administration; Writing—Review & editing. Andy R. McKinley: Conceptualization; Funding acquisition; Resources; Supervision; Writing—Review & editing. Charan Ranganath: Conceptualization; Funding acquisition; Resources; Supervision; Writing—Original draft; Writing—Review & editing.

This project was supported by the Air Force Research Laboratory (AFRL) under the Human Interface and Research Technology program (https://dx.doi.org/10.13039/100000005), grant number: FA8650-14-D-6500; and by a Vannevar Bush Faculty Fellowship (ONR grant N00014-15-1-0033) to C.R. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the U.S. Department of Defense.

A retrospective analysis of the citations in every article published in this journal from 2010 to 2020 has revealed a persistent pattern of gender imbalance: Although the proportions of authorship teams (categorized by estimated gender identification of first author/last author) publishing in the Journal of Cognitive Neuroscience (JoCN) during this period were M(an)/M = .408, W(oman)/M = .335, M/W = .108, and W/W = .149, the comparable proportions for the articles that these authorship teams cited were M/M = .579, W/M = .243, M/W = .102, and W/W = .076 (Fulvio et al., JoCN, 33:1, pp. 3–7). Consequently, JoCN encourages all authors to consider gender balance explicitly when selecting which articles to cite and gives them the opportunity to report their article's gender citation balance. The authors of this article report its proportions of citations by gender category to be as follows: M/M = .523, W/M = .233, M/W = .116, and W/W = .128.