Network-level mechanisms underlying effects of transcranial direct current stimulation (tDCS) on visuomotor learning.

Transcranial direct current stimulation (tDCS) is a non-invasive brain stimulation approach in which low level currents are administered over the scalp to influence underlying brain function. Prevailing theories of tDCS focus on modulation of excitation-inhibition balance at the local stimulation location. However, network level effects are reported as well, and appear to depend upon differential underlying mechanisms. Here, we evaluated potential network-level effects of tDCS during the Serial Reaction Time Task (SRTT) using convergent EEG- and fMRI-based connectivity approaches. Motor learning manifested as a significant (p<.0001) shift from slow to fast responses and corresponded to a significant increase in beta-coherence (p<.0001) and fMRI connectivity (p<.01) particularly within the visual-motor pathway. Differential patterns of tDCS effect were observed within different parametric task versions, consistent with network models. Overall, these findings demonstrate objective physiological effects of tDCS at the network level that result in effective behavioral modulation when tDCS parameters are matched to network-level requirements of the underlying task.

Introduction

Transcranial Direct Current Stimulation (tDCS) is a non-invasive brain stimulation approach in which low level currents are applied across the scalp to influence underlying brain function (rev. in (Brunoni et al., 2014)). Although potential effects were known from classical studies, a recent reemergence of tDCS research began ~20 years ago (Nitsche and Paulus, 2000) with the demonstration that low-level currents applied over the scalp effectively modulated underlying cortical excitability as measured using transcranial magnetic stimulation (TMS) (Nitsche and Paulus, 2000). As a result, tDCS, along with other non-invasive brain stimulation approaches, have high potential to enhance cortical function (Nissim et al., 2019), enhance neurorehabilitation (Sánchez-Kuhn et al., 2017), and reverse deficits in neuropsychiatric disorders (Szymkowicz et al., 2016) such as depression (Bennabi and Haffen, 2018) or schizophrenia (Brunelin et al., 2012). Nevertheless, despite an extensive literature on tDCS effects, significant issues related to its mechanism of action remain unanswered (for discussion see Polanía et al. (2018)). In particular, models of tDCS have focused extensively on its effects on local excitability within cortex directly underlying the stimulation electrodes (Polanía et al., 2018; Bikson et al., 2018; Jacobson et al., 2012). By contrast, more recent studies focus more on network-level effects across cortical regions (To et al., 2018). Here, we used multimodal imaging combined with a well-studied (rev in Buch et al. (2017)) motor learning paradigm -the serial reaction time task (SRTT) - to investigate mechanisms underlying tDCS effects.

In the SRTT, subjects make a series of button presses based upon visual cues (Fig. 1A). When sequences are unpredictable (i.e. “random” sequences), reaction time (RT) remains constant over the course of repeated trials. By contrast, when the sequence repeats predictably (“fixed” sequences), individuals typically show a progressive reduction in RT over repeat trials even if they are not told of the sequence in advance, a process termed implicit motor learning. Motor learning in this task is distributed across a range of levels in the brain, beginning with perceptual and posterior parietal nodes of the dorsal visual system (Ashe et al., 2006) and projecting to motor cortex and premotor cortex/supplementary motor area (PMC/SMA) (Hardwick et al., 2013 b; Focke et al., 2017). Sequence information may be specifically encoded within dorsal stream visual regions (Hardwick et al., 2013b; Focke et al., 2017; Kantak et al., 2012; Savic and Meier, 2016; Grafton et al., 1998; Keele et al., 2003) as well as visual sensory regions (Gavornik and Bear, 2014b, 2014a).

Fig. 1.

Task, tDCS montages. A. Schematic Task description: Top Right: During the Serial Reaction Time Task (SRTT) subjects were cued to press on each trial on a key corresponding to 1 of 4 locations on the display screen. Crosses, corresponding to each of the four designated keys on the keyboard, were persistently present on the monitor. Subjects also performed a Simple Reaction Time (SRT) task, in which a stimulus was presented with jittered ISI at a single location and only a single response type was required. Left: The figures present the sequence of task presentations and their relationship to tDCS stimulation conditions. For EEG all four tDCS conditions were implemented on separate days. For MRI only sham and visual-cathodal were implemented. The 30 s periods indicate the ramp up and down of the tDCS stimulation window. B. tDCS montage: Pads were placed as shown for the Motor-cathodal/anodal (top) and Visual-cathodal (bottom) stimulation conditions, leading to electric field (EF) intensity distributions affecting primarily anterior and posterior brain regions, respectively. As the EF intensity maps do not show direction of current flow the maps for the anodal and cathodal motor stimulation conditions would be the same.

In general, task-dependent modulation (Boroda et al., 2020) of motor activity within this extended visuomotor networks, including in the SRTT is reflected in alterations in coherence within the β (14–24 Hz) frequency range (Brovelli et al., 2004; Schoffelen et al., 2008; Hardwick et al., 2013a), as well as in fMRI functional connectivity between regions (Tzvi et al., 2014). To our knowledge, however, a multimodal approach utilizing these measures has not previously been used to assess tDCS effects on motor learning, in part due to technical limitations in the use of scalp-based β-activity to assess intracranial effects, and in part due to limitations of current methods by which SRTT behavioral data are typically analyzed. Here, we used two innovations to permit the application of the β-coherence approach to the study of tDCS effects on long-range cortical synchronization. First, instead of using mean RT across groups of trials, as is typical, we resolved RT on a trial by trial basis. Second, we used source-space β-coherence analyses complemented with fMRI connectivity measures to overcome the limitations of standard, surface-based approaches.

Traditionally, RT changes in the SRTT have been studied using mean RT values averaged across blocks of trials (Huang et al., 2018b), often with the underlying assumption that the large decrease in RT (typically > 100 ms) reflects increased efficiency of motor cortex. However, a more parsimonious explanation of the learning effect is that as a consequence of learning the sequence, individuals are increasingly able to predict where the next stimulus will appear and thus are able to prepare in advance where to press, so that the stimulus only needs to indicate when to press. By contrast, in other trials, the subject needs to wait for the stimulus and react when it occurs. In classical terms (Donders, 1969), a task in which the stimulus instructs both type of response and timing is termed a “choice reaction time” (CRT) task, whereas one in which the stimulus instructs only response timing is termed a “simple reaction time” (SRT) task. RT is longer in a CRT vs. SRT task by several hundred milliseconds, reflecting the increased processing required to first determine response type and only then to initiate the response (Donders, 1969; Greenhouse et al., 2015).

Here, we hypothesized that in the random condition, in which it was (by definition) impossible to predict where the next stimulus would appear, a single distribution of “reactive” long-RT trials would be observed. By contrast in the fixed condition, a mixture of fast, “predictive” and slow, “reactive” trials would be observed, with the reduction in mean RT reflecting primarily an increase in percentage of predictive vs. reactive responses over time. Given the large difference in RT between slow, reactive and fast, predictive responses, we further hypothesized that the different response types would depend upon differential connectivity patterns across networks, rather than alteration in response magnitude within motor cortex itself. Finally, we hypothesized that differences would be most apparent within the β-frequency range.

To test this hypothesis, we analyzed RT changes across subjects at both the mean level across trials as usual (e.g. (Buch et al., 2017; Hardwick et al., 2013b; Savic and Meier, 2016)), and the single trial level. In parallel, we collected multimodal event-related potential (ERP) and fMRI measures simultaneously to tDCS. ERP measures have high temporal resolution and thus can be analyzed on a trial-by-trial basis, permitting investigation of neural networks engaged during slow (CRT-like) vs. fast (SRT-like) responses. However, the spatial resolution is limited, so that source analyses must be confirmed by other approaches. fMRI possesses high spatial resolution, permitting precise localization of structures involved in the task. However, temporal resolution is low, preventing differentiation of activation patterns relating to different response types (i.e. fast vs. slow). By using convergent measures, we were able to obtain high temporal as well as high spatial resolution.

In ERP, cortical control over motor activity is indexed most strongly by event-related desynchronization (ERD) of ongoing beta-frequency (10–24 Hz) rhythms (Pfurtscheller and Lopes da Silva, 1999; Roelfsema et al., 1997; Gompf et al., 2017; Gladwin et al., 2008; Jasper and Penfield, 1949; Neuper and Pfurtscheller, 2001). These rhythms reflect interactions both within motor cortex and between motor cortex and related regions (rev. in (Khanna and Carmena, 2015; Weinrich et al., 2017)). Here, we used a source-space analysis approach (Scherg et al., 2019; Gross et al., 2001; Mehrkanoon et al., 2014) to evaluate pairwise β-coherence (Singer and Gray, 1995; Bressler, 1995 b; Buzsaki and Draguhn, 2004) between dorsal-visual, PMC/SMA and motor nodes of the canonical SRTT circuit (Kantak et al., 2012; Keele et al., 2003) as a function both of motor learning and tDCS effect, and used fMRI to verify the source location used for ERP analyses. In general, increased β-ERD is associated with increased BOLD response within motor cortex, suggesting that it may index bringing motor regions “online” during the task. Nevertheless, non-linear relationships between β activity and BOLD are also observed. Thus, depending on the task, both increases and decreases in β-ERD and both motor activity and BOLD response may be observed (Gompf et al., 2017).

From the fMRI data, we also evaluated effects of tDCS on effective connectivity between frontal, motor and visual regions using a general psychophysiological interaction (gPPI) approach (McLaren et al., 2012). In general, PMC/SMA neurons are involved in “decision making” whereas motor cortex neurons are more proximate to the specific response (Huang et al., 2018 b; Lara et al., 2018; Tanji and Mushiake, 1996; Nachev et al., 2008). We hypothesized that during motor learning, slow and fast responses would be differentiable based upon the relative degree to which they depend upon Dorsal-visual → PMC/SMA vs. Dorsal-visual → Motor cortex interactions. These network measures were then used to evaluate the mechanism of action of tDCS.

As in other tDCS studies of SRTT (Manoach et al., 2004; Nitsche et al., 2003), we applied both anodal and cathodal stimulation over motor cortex to evaluate polarity-specific effects on behavior and physiological activity. In addition, given our network hypothesis, we also stimulated over dorsal stream visual cortex (Pobric et al., 2018; Zito et al., 2015; Antal et al., 2004 a, 2004 b; Heinen et al., 2016), which is known to be implicated in visuomotor learning (Keele et al., 2003). These studies thus test the hypothesis that beneficial effects of tDCS are mediated at the network level, and that effects on motor speed may be obtained even without direct targeting of motor cortex.

Methods

Participants

Seventeen (4 female) neurologically normal paid volunteers, aged 19.4–54.0 (mean = 35.2, SD ± 10.8), participated in the ERP study. Eight (2 female) also participated in the fMRI section of the study in a counterbalanced order. All subjects provided written informed consent, and the procedures were approved by the Nathan Kline Institute/Rockland Psychiatric Center Institutional Review Board and ethics committee. All subjects reported normal or corrected-to-normal vision. All were right-handed.

Stimuli and experimental design

SRTT

Subjects were instructed to press one of four visually cued color-coded keys on a standard computer keyboard with the fingers of their right hand as quickly and accurately as possible following presentation of a corresponding visual cue consisting of a color-coded square corresponding to the keys (Fig. 1A). A 3-s response window was permitted following each cue. Crosses, corresponding to each of the four designated keys on the keyboard, were persistently present on the monitor. Stimuli were presented on a computer monitor (Iiyama Vision Master Pro 502, model # A102GT) located 143 cm from the subject and consisted of four crosses that collectively subtended ± 1.4 ° visual angle from the center of the screen.

The studies were performed using a 5-element repeat sequence, modeled after previous studies (Manoach et al., 2004; Nitsche et al., 2003). We used four different SRTT sequences, one for each tDCS condition (see Section 2.3), i.e. Sham (3, 1, 4, 2, 4), Motor-cathodal (2, 3, 1, 2, 4), Visual-cathodal (1, 3, 4, 2, 3), Motor-anodal (4, 2, 1, 3, 2), to ensure no one received the same sequence twice. Two blocks of SRTT, 12 min each, were administered during tDCS/EEG (Fig. 1A). Each block consisted of 12 self-paced runs, with random runs at positions 1 and 10 of the sequence, as per (Nitsche et al., 2003). A single block was repeated 10-min post-tDCS.

Motor learning was defined as the difference in RT in each fixed block minus RT in the random block. Similar methods were used in the fMRI experiment (Fig. 1A), except that SRTT was delivered through MR-compatible liquid crystal display goggles (Resonance Technology Inc., Northridge, CA).

SRT

A control task was also administered prior to, during and following tDCS. In this condition only a single (central) stimulus was presented in each trial and subjects responded by pressing a single button with their right index finger as quickly as possible following cue stimulus presentation. Stimuli were presented with a mean ISI of 100 ms. A jitter of ± 50-ms was used to prevent self-paced tapping at a fixed ISI.

Transcranial direct current stimulation

tDCS application

tDCS was applied by a saline-soaked pair of surface sponge pads (3 × 3 cm) using the battery-driven, NeuroConn DC-Stimulator MR (NeuroConn, Ilmenau, Germany). During the ERP section of the study, the participants received four stimulation conditions (Sham, Motor-cathodal, Visual-cathodal, Motor-anodal) using a constant current of 2-mA intensity applied for 30 min during the task performance (Fig. 1A). Each stimulation condition was administered on a separate day (at least 36 h apart) for each subject in counterbalanced order.

Pad placements for the Motor-cathodal and Motor-anodal conditions followed the M1-SO (left primary motor-right supraorbital) scalp positions used in prior tDCS SRTT studies (rev in Buch et al. (2017)) (Fig. 1B). Ensuring that the sponges are not so wet as to drip we applied a thin layer of Ten20 paste (Weaver and Company) to provide good contact with the scalp and prevent the sponges from moving. We then carefully placed one pad under the EEG cap at position C5 and the other at FP2 in accordance with the 10–20 EEG electrode system. We took care that no prior local skin inflammations were present at pad locations. For Visual cortex stimulation, the anode pad was placed over the vertex (Cz) and the cathode pad was placed on the scalp area (POz) overlaying the Glasser atlas-defined dorsal visual area (Glasser et al., 2016). For sham stimulation, the pads were placed in the same positions as for motor stimulation; however, the stimulator only delivered 30 s of ramp up and down. During simultaneous fMRI-tDCS experiments, only visual-cathodal and sham stimulation conditions were used. For placement of stimulation electrodes, we followed the above procedure using the damp sponges along with a layer of Ten20 paste and using the EEG cap cloth, without the EEG electrodes and holders, to hold the electrodes in place. The stimulation cables were run through the side-openings of the head-coil as recommended by the manufacturer (neuroCare Group 2016). The stimulation pad wires were each equipped with 5.6 kΩ resistors to avoid temperature increases that could result from induction voltages. For the MRI-tDCS setup we followed the precautions detailed in (Esmaeilpour et al., 2019) and closely replicated the set up in (Antal et al., 2011b). The stimulator device was connected to the subject via a waveguide with a length:width ratio of 4:1 with RF filters on either side. The RF filters, one placed in the control room at the opening of the waveguide and one placed at the head end opening of the magnet bore, proximal and off-centered from the head coil (neuroCare Group 2016), had an approximate attenuation of 60 dB within a frequency range of 20–200 MHz. The cable ran at the edges of the scan room as far as possible from the magnet bore and was weighed down using sandbags to prevent movement and curling of the cable. Impedance was checked continuously and was kept under 10 kΩ during the stimulation. The MRI machine triggered the onset of the task on the task delivery PC using a BNC connection between the two. One minute prior to the onset of the task the tDCS stimulation was started manually to allow for the ramp up of the current.

Finite-element modeling of average electric field strength was performed on the MNI-152 head (6th generation, non-linear - T1-weighted), using the ROAST (Huang et al., 2018 a) toolbox in MATLAB. Thirty empty slices were first added to the MNI volume in all six directions to allow room for simulated electrodes. Then both montages were simulated with 3 cm square pad electrodes of 3 mm thickness, using standard conductivity values for pads and gel. The motor stimulation montage comprised pad electrodes centered at C5 and Fp2, passing − 2 mA and + 2 mA, respectively (motor-cathodal in this case). The visual stimulation montage comprised electrodes centered at POz and Cz, passing −2 mA and + 2 mA, respectively. Electrical field maps output by ROAST as NIfTI volumes were then mapped onto a standard HCP averaged surface (WU-Minn HCP Data - 500 Subjects dataset, 32k mesh resolution) in HCP Workbench (Fig. 1B), and field strength averages were computed for each ROI in the Glasser parcellation (Glasser et al., 2016).

tDCS discomfort

We used the Faces Pain Scale (Hockenberry et al., 2005) to measure discomfort caused by tDCS application after each session. This scale consists of 6 faces with face 0 indicating a happy face denoting “No Hurt” and face 5 indicating a crying face denoting “Hurts Worst”. This scale illustrates physical pain and is easy to understand. We also asked the participants to verbalize what type of hurt they felt if any i.e. itching, burning, pain, etc.

EEG and fMRI data acquisition

EEG acquisition

Continuous EEG was acquired through Brainvision Brainamp MR Plus amplifier system using 32 scalp active electrodes, impedances < 5 kΩ, referenced to the FCz electrode, bandpass filtered from 0.05 to 100 Hz, and digitized at 500 Hz. Data were re-referenced to average-reference and analyzed offline using BESA Research, version 5.3 (Brain Electric Source Analysis, BESA GmbH), EEGLAB (Delorme and Makeig, 2004), ERPLAB (Lopez-Calderon and Luck, 2014) and Matlab software, version 2017a (MathWorks). The recording of EEG during tDCS requires particular attention because neurostimulation may also affect electrophysiological measurement (Gebodh et al., 2019). The EEG system used here, having DC amplifiers, eliminates the delayed response after amplifier saturation that is observed with EEG systems relying on AC recordings, which contribute to data clipping and recording artifacts. The system provides a dynamic range of 110 dB that handles the offset potentials during tDCS application without the need of a high-pass filter and provides equal frequency responses at all channels, aiding to mitigate data distortion during stimulation (Levit-Binnun et al., 2010).

Prior studies have successfully removed EEG artifacts caused by concurrent tDCS (Roy et al., 2014). These include cardiac and ocular motor distortion as inherent and motion and myogenic distortion as noninherent physiological artifacts (Gebodh et al., 2019). The non-inherent physiological artifacts are greatly reduced by instructing the participants to minimize their movements and to relax their neck and jaw muscles. Nevertheless, tDCS can accentuate these artifacts if present.

Here we applied an artifact criterion of ± 70 μV at all scalp sites, in addition to manual inspection, to reject epochs contaminated by these artifacts. The inherent physiological cardiac artifact is not event-related and can be baseline corrected. As previously suggested (Gebodh et al., 2019) baseline correction might not be adequate to eliminate this artifact. In the frequency domain however, the peak cardiac artifact signature is within 1–1.2 Hz frequency range which is below our frequency range of interest i.e. 10Hz-24 Hz (see results). For detection and removal of ocular motor and other artifacts presenting with characteristic topographic signatures we employed blind source separation using Generalized Eigenvalue Decomposition as spatial filtering approach (Cohen, 2017; Comon, 1994; Parra and Sajda, 2003; Parra et al., 2005). The average EEG epoch acceptance rate was 85%.

fMRI acquisition

Scanning was conducted on a 3T Siemens TiM Trio system (Erlangen, Germany), using a standard twelve-channel circularly polarized head coil, at the Center for Biomedical Imaging and Neuromodulation at NKI using standard measures (Sehatpour et al., 2010; Sehatpour et al., 2006) and considerations for concurrent MRI-tDCS implementation outlined in (Esmaeilpour et al., 2019). Briefly, we acquired two resting state scans, one during sham and one during active tDCS, in axial orientation aligned to the AC-PC plane (TR = 2500 ms, TE = 30 ms, FOV = 216 mm, matrix = 72 × 72, 38 3.3 mm slices, 0.7 mm gap, acceleration factor = 2) with a time series duration of 240 TRs and task-based (same sequence) scans with a time series duration of 209 TRs during task + tDCS stimulation. Blood oxygenation level dependent (BOLD) data were acquired during resting state and task. Motor responses were registered using a custom-made MR-compatible keyboard. T1-weighted structural scans also were collected (TR = 1900 ms, TE = 2.52 ms, TI = 900 ms, FOV = 256 mm, matrix = 246 × 256, 176 1-mm slices, no gap, acceleration factor = 2) aligned to the AC-PC plane.

Data analysis

Behavioral data analysis

Mean RT analysis:

For initial analyses, RT data from fixed blocks were log-transformed, averaged across trials within a block, and normalized by block relative to mean RTs from the interspersed random blocks. Exponential decay functions were calculated for each condition using group mean data using PRISM Graphpad 7.0 (GraphPad Software, San Diego, CA). Values of the plateau and time-constant values were tested across using sum-of-squares F-tests across the four stimulation conditions. If a significant difference was observed across conditions, pairwise comparisons between active tDCS conditions and sham were performed. Confirmatory analyses were conducted using mixed-model regression with condition as a factor and run number as a covariate. Post-hoc comparisons were conducted using Sidak correction for multiple comparisons. Potential order effects were assessed using session number for each condition as an additional factor.

Single-trial RT distribution analyses:

For follow-up analyses, single-trial log-RT distributions were compared across conditions using single vs. dual-Gaussian models using GraphPad 7.0 non-linear curve fitting functions. For each analysis, both single Guassian and dual Gaussian fits were considered. For the single Guassian model, data were fit to the formula

where Total number of responses, Mean, and SD were modeled parameters. All parameters were constrained to be positive values. For the dual Gaussian model, a second set of parameters corresponding to the second Gaussian distribution were added. Starting values were provided based upon apparent peaks in the histogram plots.

In all cases, the simpler model (single Gaussian) was chosen unless the more complex model (Dual Gaussian) was shown to be statistically superior. Comparison between models was assessed using a goodness-of-fit ANOVA. In addition, absolute goodness of fit (R2) was required to be > 95% for all accepted fits. Initial analyses were performed using data from the random repeat runs only, which yielded unimodal models in all cases. Mean RT values from the random runs were used to constrain the slow RT component for the subsequent analyses of RT data from the fixed-sequence blocks. Initial values for each model were provided based upon visual inspection of RT histograms. Analyses were conducted both by quarter to evaluate stability of RTs over the course of the training and collapsed across quarters to compare coefficients.

Comparison of percentage fast responses across tDCS conditions was performed by comparing dual Gaussian models in which the ratio between fast and slow responses was assumed to be constant vs. those in which it was assumed to vary across conditions. The simpler model (all percentages equal) was accepted unless the more complex model (percentages different) was found to be statistically superior. For individual runs, because there were not sufficient trials to model single- vs. dual-Gaussian fits, a cut value of log RT = 2.35 log-ms (223.9 ms) was chosen to best discriminate fast from slow responses. For display purposes, the absolute number of trials across all RTs was normalized to 100% for all conditions.

EEG data analysis

Sensor-level analysis:

Epochs of continuous scalp-recorded EEG extending from 200 ms before to 200 ms after each motor response were used to compute the response-locked ERP. The data from the two missing electrodes due to pad placements in each stimulation condition was replaced using Spherical Spline Interpolation (SSI) (Perrin et al., 1987). Baseline was defined as the 200 ms to 100 ms pre-response interval. To obtain the envelope amplitude and the phase of a specified frequency band as a function of time, time-frequency transformation was performed by using complex demodulation (Hoechstetter et al., 2004; Sehatpour et al., 2008) for frequencies of 4–50 Hz in the time-window between −500 to 600 ms. Baseline was defined as the 200 ms to 100 ms pre-response interval. Frequencies were sampled in 2-Hz steps; latencies were sampled in steps of 25 ms. This corresponds to a time-frequency resolution of +/− 2.83 Hz and +/− 39.4 ms at each time-frequency bin (full width at half maximum).

β-ERD values were calculated using temporal spectral evolution (TSE) defined as the relative power change at a time-frequency bin compared with the mean power over the baseline epoch for that frequency (Pfurtscheller and Lopes da Silva, 1999; Sehatpour et al., 2008). A one-way ANOVA with post-hoc Sidak correction was performed to test for a significant differences of TSE values between sham and each active condition at the C3 scalp electrode site overlaying the left motor cortex. The corrected significance level α was set to 0.05.

Source-space analysis:

Intracranial sources of beta-activity were assessed using a Beamformer approach (Van Veen et al., 1997; Sekihara et al., 2001), which involves the following steps: 1) For each channel single-trial data in time domain is transformed into time-frequency domain in order to compute the complex time-frequency signal (Sehatpour et al., 2008); 2) Complex cross-spectral density matrices is then computed for each trial. 3) A forward model is applied, and a lead-field matrix is estimated; and 4) The brain in Talairach space is divided into a grid and the normalized task-dependent contribution (q) to a given time-frequency range of interest from every location on this grid is then estimated.

Each channel single-trial data in time domain is transformed into time-frequency domain, deriving as the complex time-frequency signal of channel x in trial n at frequency f and latency t, characterized by its amplitude A and its phase φ.

Prior to computation of the Beamformer we removed the event-related signal from each trial for calculations of the estimates in the frequency domain (Sehatpour et al., 2008; Van der Lubbe et al., 2016). This transformation ensures stationarity in the EEG data notably during concurrent tDCS application.

For each trial n the regression coefficient b n is computed as: where

From each trial n, the correlated fraction bn of the time-frequency transform of the average waveform was subtracted. This approach accounts for the amplitude fluctuations between trials and ensures stationarity.

Complex cross spectral density matrices C are then computed for each trial: Here, * indicates the complex conjugate.

The output power P of the beamformer for a specific brain region at location r is then computed as (Gross et al., 2001): Here is the inverse of the average of C (f, t) over trials and the time-frequency range of interest, L is the leadfield matrix (i.e. the magnitude of the signal each source contributes to each recording sensor (Scherg and von Cramon, 1985) of the model containing the regional source (Sehatpour et al., 2006; Scherg and Picton, 1991) at location r. In computing the lead-field matrix we used the standardized finite element model (FEM) implemented in BESA which has been created from an averaged head using 50 individual MRIs in Talairach space. The FEM model provides a realistic approximation to the averaged head and uses three compartments: brain/CSF, skull and scalp to describe the electrical conductivity distribution inside of the head. T is the matrix transpose and tr’ [] is the trace of the [3 × 3] submatrix (of the expression in the bracket) for the source at location r.

We then normalized the power P(r) with the power at the corresponding time-frequency interval of the sham condition Psham (r) to obtain the value q(r):

The brain in Talairach space is divided into a grid with a resolution of 5 mm3 and the beamformer image is constructed from values q(r) computed from every location on this grid. q values are then shown in% where q[%] = q * 100. This image is then extrapolated to a resolution of 1 mm3 and projected to an inflated brain image derived from an MRI of equal resolution. Since in the computation of beamformer image regional sources having three orthogonal vectors (i.e. radial, tangential and oblique) are used, projection onto an inflated brain surface more accurately represents the spread of the cortical activation. The cortical regions with the highest q value are then seeded with a virtual source (Scherg et al., 2019; Hoechstetter et al., 2004), revealing three distinct cortical regions at MNI coordinates: PMC/SMA [−36, 2, 59] (left BA-6), Motor [−43, −20, 60] (left BA-4), Visual [−38, − 87, 9] (left BA-18). The goodness-of-fit value of the model (% variance explained) is 90.0 SE ± 2.0. Seeded sources allow for individualized fitting of the dipolar angle accounting for variations in subjects’ cortical anatomy (Scherg et al., 2019). Source mean β-frequency ERD amplitudes for each identified generator (virtual source) were determined on a single-trial basis across all tDCS conditions and evaluated by ANOVA with factors of tDCS condition and region.

In the next step single-trial source TF data was used to derive coherence measures across this cortical network as a measure of functional connectivity (Bressler, 1995 a; Bressler et al., 1993; Fries, 2005).

Coherence ranges from 0 (no coherence) to 1 (maximum coherence). To determine the probability that coherence at a particular timefrequency sampling point is significantly higher than what is expected from random fluctuations is investigated based on an approach suggested by Lachaux et al. (1999) and previously implemented and described by our group (Sehatpour et al., 2008).

To investigate the probability that the coherence in sham differed significantly from the coherence in each of the other conditions, and to determine if two active conditions differ, the individual subject mean coherence estimates were then subjected to a permutation cluster analysis (Bullmore et al., 1999; Maris and Oostenveld, 2007). This approach is carried out in two general steps.

In the first step a Student’s paired t-test is carried out for every time-frequency (TF) bin to determine if there is a significant difference between the two conditions in the group. Here a cluster alpha level of 0.05 is set which allows us to identify the TF bins whose t-values exceed the 95th quantile threshold and can be considered as a candidate to be included in a cluster of TF bins, which in our case are based on temporal and spectral contiguity.

A cluster value is obtained by summing the t-values of the individual data bins in a cluster. This value serves as a test statistic for the next step of the analysis. In addition to the results obtained in this preliminary parametric step of the analysis we considered clusters to be submitted to the next step of the analysis, for further interrogation by permutation testing, based on our a-priori hypothesis (see the introduction) as well as the results observed from the scalp data which indicated the modulatory effects of tDCS within a time window of 100 ms prior to the motor response in the β frequency range 10Hz − 24HZ.

In the second step of the analysis the clusters within the a priori TF window are subjected to permutation testing wherein the coherence data for sham gets systematically interchanged with the coherence data of the test condition. For each permutation, a new t-test is obtained per TF bin and a new test statistic (cluster-level summed t-values) is computed.

Here we have used 2000 permutations (drawn randomly without repetitions) from all possible permutations, i.e., 217. From the distribution of the test statistics obtained from our permutations we then calculate the proportion of the test values that are larger than the value obtained from the initial cluster obtained in step 1. Hence if less than 5% of all values are larger than the initial test value it is assumed that the data of the two conditions are not interchangeable with a chance level greater than 95% i.e. (P <.05).

fMRI data analysis

Activation analysis (BOLD).

fMRI data processing was carried out using FEAT (FMRI Expert Analysis Tool) Version 6.00, part of FSL (FM-RIB’s Software Library, www.fmrib.ox.ac.uk/fsl). The following pre-statistics processing was applied following the removal of the first three TRs; motion correction using MCFLIRT (Jenkinson et al., 2002); non-brain removal using BET (Smith, 2002); spatial smoothing using a Gaussian kernel of FWHM 6.0 mm; grand-mean intensity normalization of the entire 4D dataset by a single multiplicative factor; high pass temporal filtering (Gaussian-weighted least-squares straight line fitting, with sigma = 45.0 s). Time-series statistical analysis was carried out using FILM with local autocorrelation correction (Woolrich et al., 2001). Registration to high-resolution structural and standard space images (2 mm MNI space) was carried out using FLIRT (Jenkinson et al., 2002; Jenkinson and Smith, 2001).

To examine effects of stimulation, first level scans for the Sham condition were entered into a mixed-effects (OLS) analysis. To examine condition effects, first-level scans for Sham and Visual conditions were entered into a repeated-measures analysis. Contrasts of interest were the sham-vertex and vertex-sham comparisons. Z (Gaussianized T/F) statistic images were thresholded using clusters determined by Z > 2.3 and a corrected cluster significance threshold of P <.05, with minimum cluster size determined according to Gaussian Random Fields (GRF) theory (Worsley, 2001).

RSFC:

Functional connectivity (rsFC) was computed from resting-state fMRI scans during sham and active tDCS, using DPARSFA as described previously (Hoptman et al., 2014; Chao-Gan and Yu-Feng, 2010). The first 10 images for each time series were discarded. Next, images from each session were motion-corrected. The T1-weighted image was registered to the functional data, linear trends were removed from the motion-corrected data, nuisance covariates were applied (24 motion parameters, global, CSF, and WM signals), registered to the EPI template in 3 mm MNI standard space, and smoothed with a 6-mm FWHM kernel. Images were then filtered using a 0.01–0.1 Hz bandpass.

For assessment of sham/active (visual-cathodal) rsFC changes, a spherical region of interest (ROI; 5 mm radius; MNI coordinate −42, −18, 54, left pericentral gyrus) was created from the peak location of motor activation on the sham activation images across subjects. This was used as a seed for the rsFC analyses. rsFC values were converted to Z scores using Fisher’s r-to-z transform. These data were analyzed using repeated-measures mixed-effects (OLS) analysis. Statistical images were thresholded using clusters determined by Z > 2.3 and a (corrected) cluster significance threshold of P=.05 (Worsley, 2001).

gPPI analysis.

A generalized context-dependent psychophysiological interaction analysis (gPPI) (McLaren et al., 2012; Cisler et al., 2014) was implemented to examine task-associated changes in functional connectivity in relation to tDCS stimulation. Analyses were carried out with AFNI software (Analysis of Functional NeuroImages, http://afni.nimh.nih.gov/afni/) following similar pre-processing steps as described above for RSFC. Seed time series from PMC/SMA, motor, dorsal and ventral regions were extracted from the pre-processed fMRI data during sham and cathodal visual tDCS conditions. Seed regions were based on the location of peak activation during sham fMRI across all subjects. Each extracted time series was detrended and deconvolved with the gamma variate hemodynamic response function (HRF).

For each seed, the PPI regressor for the task was generated as the product of the deconvolved-seed time series and a regressor based on the timing of the task. Separate individual-level gPPI general linear model (GLM) analyses were then conducted for each seed using AFNI 3dDeconvolve. Each model included the PPI regressor, the seed time series and the regressor of the original model. Including both seed and task regressors of the original model in the GLM allowed us to test the PPI term as the interaction effect above and beyond the main effects. The resulting beta coefficient for the PPI regressor is thus an estimate of the amount of signal explained by both the response in the seed region and the stimulus condition. Finally, mean PPI beta coefficients (for each seed and tDCS condition) were extracted from each seed region in a pairwise fashion and entered into group-level analyses.

Results

tDCS results

Field strength mapping

As expected, stimulation over motor cortex produced highest field strength within left premotor (0.320 ±.038 V/m) and somatomotor (0.315 ±.038 V/m) regions, with lower mean field strength within visual regions (mean 0.102 ±.008 V/m). By contrast, visual cortical stimulation produced greatest inward current flow within the dorsal visual region (0.219 ±.028 V/m) and outward current flow within the superior parietal region (0.219 ±.031 V/m) (Fig. 1B).

Discomfort

All subjects tolerated the tDCS/EEG and tDCS/MRI well. None of the participants reported a score higher than 2 on the scale of 0–5 (during EEG: mean = 1.2, SD ± 0.6), (during MRI: mean = 1.1, SD ± 0.4). No correlation between discomfort level and tDCS condition (including sham) was found (r = 0.036; p=.892). The subjective feeling of discomfort was primarily reported as itching (64%), tingling (24%) and burning (12%). No adverse events occurred during this study.

SRTT results

SRTT behavior

The effects of tDCS on motor learning were assessed using exponential modeling of behavioral data across runs.

Fixed condition:

During stimulation, RT declined according to a single exponential across runs, (Fig. 2A, bottom tracings) with a significant difference in the plateau across conditions (F3,72 = 15.91, p <.0001). Motor cathodal (F1,36 = 10.8, p=.0023) and Visual (F1,36 = 26.93, p <.0001) stimulation both induced significant RT differences vs. sham, whereas Motor-anodal stimulation was statistically ineffective (F1,36 = 2.13, p=.15) (Fig. 2B, top). Dorsal visual stimulation also increased the rate of improvement (F1,36 = 7.62, p=.009), whereas other conditions did not affect motor learning rate (both p >.6). In the random condition there was no significant change over time in any group (all p >.4), or significant effect of tDCS (F3,12 = 1.736, p=.21) (Fig. 2A, top tracings).

Fig. 2.

Behavioral effects in the SRTT during tDCS. A. Mean RT by tDCS condition across fixed and random runs both during tDCS stimulation and 15-min post-stimulation. RTs during random runs are shown in the top part of the graph and remain constant over runs. RTs during the fixed condition are shown in the bottom part and show exponential reduction over time. B. Mean final RT by tDCS condition during (top) and after (bottom) stimulation. Values are the plateau values from panel A. C. RT distribution by quarters for Random (top) and Fixed (bottom) sequences, showing that responses fit to unimodal distributions across for the random condition across all quarters vs. bimodal distributions in the fixed condition. Q1 represents combined responses for runs 1–6, Q2: runs 7–12; Q3: runs 13–18; and Q4: runs 19–24. Thus, in each quarter there is one random run and 5 fixed runs. Responses from after stimulation are combined across all runs (1–12). Note log RT axis. D. RT distribution collapsed across quarters during the fixed sequence by tDCS condition showing bimodal fits under all conditions. RT distributions for the random sequences (not shown) were unimodal under all tDCS conditions. E. Percent fast responses by run: Fast responses were defined as those with RT < 2.35 log-ms (223.9 ms), showing progressive increase with run number in the fixed condition. The number of fast responses during the random runs (not shown) was not statistically different from 0 for any condition. Inset: Plateau value for%fast responses across conditions. * * p <.01 vs. sham; * * * p <.001.

Similar results were obtained using a mixed-model regression analysis with factors of stimulation condition and with run number as a covariate. As in the exponential fit model, there was a highly significant effect of stimulation condition (F3,1549 = 21.3, p <.0001). Difference between sham and both Motor cathodal (p=.013) and Visual (p <.0001) stimulation remained significant following Sidak correction for multiple comparisons.

When behavior assessments were repeated 15-min after the end of stimulation (Fig. 2A, right), no significant further improvement occurred in any condition. The mean RT across blocks was also significantly different across conditions (F3,36 = 5.65, p=.003), with significant differences for Motor-cathodal (F1,18 = 13.8, p=.0016), Visual (F1,18 = 6.81, p=.02) and Motor-anodal (F1,18 = 9.53, p=.006) vs. sham (Fig. 2B, bottom). By contrast, no significant tDCS effects were observed for the random condition (all p >.4).

Single-trial RT distribution:

For assessment of single trial RT distribution, data were first log transformed for normalization and then split into quarters, corresponding to 1st half of block 1, 2nd half of block 1, 1st half of block 2, and 2nd half of block 2. Thus, each quarter contained data from 5 fixed runs and 1 random run. This permitted assessment of evolution of response profile across quarters. Subsequent analyses used data collapsed across quarters to compare across tDCS conditions.

During random runs, single unimodal distribution was observed across all 4 quarters in the sham condition (Fig. 2C), with mean RT across all quarters of 2.676 ± .002 log-ms (474.2 ± 1.0 ms). Similar unimodal distributions were observed for all quarters during active stimulation as well (Supplementary Table 1). In all cases, R2 goodness of fit to a single normal distribution was > 0.99.

By contrast, the RT distribution in the fixed condition fit preferentially to a bimodal (2-Gaussian) distribution across all 4 quarters in the sham condition (Fig. 2C). When comparisons were performed across tDCS conditions bimodal fits were statistically superior to a unimodal fit across all quarters and stimulation conditions (F18,478 = 58.3, p <.0001). Moreover, in all cases R2 goodness of fit to the bimodal distribution was > 0.95 (Supplementary Table 1). When data were collapsed across quarters, unimodal fits were again observed for all random conditions (not shown), and bimodal fits were again observed for all fixed conditions (Fig. 2D), with increased percentage fast responses during Motor-cathodal and Visual stimulation vs. sham (Supplementary Table 1).

When data were analyzed by run using a cutoff value of fast vs. slow responses of 2.35 log-ms (223.9 ms), the percent of fast responses increased exponentially across the initial blocks (Fig. 2E). As with RT, the plateau value for percentage fast responses was significantly different across conditions (F3,72 = 10.62, p <.0001), with highly significant effects of Motor-cathodal (F1,36 = 21.32, p <.0001), Visual (F1,36 = 30.5, p <.0001) and a small, but significant, effect of Motor-anodal (F1,36 = 4.62, p=.038).

Following stimulation, bimodal distributions were again observed during the fixed sequence under all stimulation conditions (all p <.002), with a mean RT for the fast component of 2.135 ±.007 log-ms (136.8 ms). The difference in percentage fast responses was statistically reliable across conditions (F3,100 = 2.95, p=.04), with an increase in percentage fast trials for both Motor-cathodal (F1,50 = 5.62, p=.02) and Visual (F1,50 = 4.16, p=.047) stimulation relative to Sham, whereas no significant difference was observed for Motor-anodal stimulation (F1,50 = 2.53, p=.12).

Order effects:

Visual stimulation was added later in the study and thus was the last session for most subjects. Nevertheless, there was no significant effect of order (F3,1444 = 1.33, p=.26) or order X condition interaction (F5,1443 = 1.22, p=.3). Furthermore, the condition effect remained significant (F3,1448 = 4.96, p=.002) with significant post-hoc difference between visual and sham stimulation (LSD p=.04).

SRTT neurophysiology

Scalp neurophysiology:

Neurophysiological responses were back averaged from the motor response. As expected, a premotor potential was observed ~25 ms prior to the response (Fig. 3A). tDCS modulated both the topography (Fig. 3A) and amplitude (Fig. 3B) in a location- and polarity-dependent fashion.

Fig. 3.

Average scalp β activity: A. Scalp topography of SRTT: response-locked motor potential at 25 ms prior to the key press showing modulation of the event-related topography under the various tES conditions. The green circle represents the C3 electrode scalp site overlaying the motor cortex. B. Response-locked ERP: observed at C3 showing tES modulation of the scalp potential shortly prior to the key-press (zero). Bottom panel: Mean amplitude of the ERP from the shaded area − 75 ms to − 25 ms. Main effect of condition: F3,89,161 = 3.595, p=.013; Motor-cathodal vs Sham: mean diff= 0.170, std err = 0.060, p=.013; no other significant differences were observed. C. Response-locked TSE: power spectrogram at C3 showing the tES modulatory effect in the time-frequency (TF) domain. The bounding box indicates TF window of 10–24 Hz at −75 to − 25 ms prior to key press (zero). D. Response-locked evoked power: averaged across frequency range 10–24 Hz. Bottom panel: Mean power from the shaded area − 75 ms to − 25 ms. Main effect of condition: F3,69,208 = 6.453, p <.0001; Motor-cathodal vs Sham: mean diff= − 0.02, std err = 0.008, p=.012; Visual vs Sham: mean diff= − 0.024, std err = 0.009, p=.005; no other significant difference was observed.

In time-frequency analyses, (Fig. 3C), the premotor potential corresponded primarily to a β-frequency response during the-75 to −25 ms premotor period. tDCS significantly (p <.05) modulated the premotor β-activity, with significant effect of both Motor-cathodal and Visual stimulation relative to sham (Fig. 3D).

Single-trial β-analyses:

For these analyses, β-ERD during the −75 to −25 ms premotor interval was projected into source space and analyzed across trials. Separate β-sources were resolved for left PMC/SMA, Motor and Visual regions (Fig. 4A). tDCS significantly modulated β-ERD amplitude at these intracranial sources in a tDCS-type-(F3,99,759 = 138.6, p <.0001) and region- (F2,99,759 = 22.9, p <.0001) dependent fashion. The Region X Condition interaction was also significant (F6,199,518 = 77.3, p <.0001). Post-hoc pairwise comparison between sham and the other tDCS conditions, using Sidak multiple comparison correction (Šidák, 1967), showed a significant (p <.001) modulation of β-ERD at each region (Supplementary Table 2).

Fig. 4.

Single-trial, source-based β analyses: A. Source-space : as obtained by Beamformer inverse solution procedure, showing locations of the PMC/SMA, motor (in convexity) and visual sources. Data are thresholded at p <.05, corrected. Inset: virtual source montage for the β-sources. B. Coherence measures under each stimulation condition: Left: (top row) for each source pair significant coherence measures (mean ± sem) within the beta-frequency range (10–24 Hz) in the 100-ms pre-response window are shown. The baseline coherence measures across the sources pairs was 0.01–0.09 (0.05 ± 0.04). These baseline measures did not differ significantly between the tDCS conditions. (bottom row) The coherence values for each time and frequency bin within this window (data cluster) was entered into a permutation test to assess the relationship between coherence levels in each stimulation condition vs. sham. * p <.05. Right: Bar graphs representing coherence measures under each tDCS condition across the intracranial source pairs. * p <.05 vs sham. C. Coherence measures slow vs fast trials: (top) Alterations in mean coherence values across fast and slow SRTT trials for each source pair under sham condition. (bottom) Slow vs Fast trials * p <.05. D. Coherence measures random vs fixed trials: significant differential coherence measures when comparing random vs fast SRTT trials and random vs slow SRTT trials indicating additional recruitment of PMC/SMA-Motor connectivity under the random condition. * p <.05, * * p <.01.

β-coherence was then used as an index of functional connectivity (FC) across each pair of the above three cortical sources as a function of task-type (fixed/random) and tDCS condition. Under sham stimulation, task-dependent coherence was observed across the three regions during the −75 to −25 ms interval. β-coherence was significantly larger during the pre-response interval than during baseline (0.05 ±.04), reflecting a functional brain network involved in the performance of SRTT (Fig. 4B, left). Furthermore, under sham condition, coherence estimates for fast vs. slow trials, revealed a significantly greater FC for visual cortex-motor cortex interaction in slow than fast trials, suggesting a greater reliance on this connection (Fig. 4C).

Analysis of variance for the effect of stimulation on the coherence measures showed (F3,64 = 25.3, p <.0001) tDCS modulated coherence in a location- and polarity-specific pattern (Fig. 4B, right). Post-hoc pairwise comparison between tDCS conditions, using Sidak multiple comparison correction showed, relative to sham, Motor-cathodal stimulation reduced coherence between PMC/SMA and visual cortex (0.13 ± 0.03, p <.0001) mimicking a pattern that is observed between slow vs fast trials (0.11 ± 0.03, p=.04) under the sham condition (Fig. 4C), consistent with its beneficial effect on motor learning. Visual-cathodal stimulation reduced coherence between motor and visual cortex (0.16 ± 0.03, p=.03) suggesting an alternate beneficial connectivity shift. The connectivity pattern under motor-anodal condition did not differ from sham suggesting the relatively smaller improvement in the RT (Fig. 2E) could be the result of local motor effects of anodal tDCS rather than network interactions.

Comparison of the coherence difference between motor-cathodal vs visual-cathodal stimulations showed significantly lower coherence between PMC/SMA and visual cortex (− 0.13 ± 0.04, p=.008) and significantly higher coherence between motor and visual cortex (0.18 ± 0.04, p=.02). Comparison between motor-anodal vs visual-cathodal stimulations showed significantly higher coherence between motor and visual cortex (0.17 ± 0.04, p=.03). Comparison between motor-anodal vs motor-cathodal stimulations showed significantly higher coherence between PMC/SMA and visual cortex (0.12 ± 0.04, p=.002).

Comparison of coherence for random vs fixed SRTT sequences under sham showed increased coherence in random trials when compared to fast trials (0.14 ± 0.03, p=.02) and slow trials (0.21 ± 0.03, p=.008) only between PMC/SMA and motor cortices (Fig. 4D).

SRTT fMRI

A subgroup (n = 8) of subjects also participated in fMRI scanning during the SRTT while receiving Sham or Visual (cathodal) tDCS (Fig. 5A). Consistent with the results of the β-source localization (Fig. 4A, Supplementary Table 3), significant activations were observed in left PMC/SMA along with adjacent frontal eye fields as well as in left motor cortex. When relative activation patterns were compared across brain regions, a highly significant correlation was observed between the β-source and fMRI measures (r = 0.81, p <.0001) (Fig. 6).

Fig. 5.

Effects of tDCS on SRTT-related fMRI: (A) fMRI BOLD activation: during SRTT vs. rest, showing significant activation in left premotor cortex/supplementary motor area (PMC/SMA), Motor cortex, and ventral/dorsal visual regions. Data were thresholded at p <.001, corrected. B. Generalized psychophysiological interaction: gPPI strength between indicated seed and target regions. C. fMRI-based resting state functional connectivity: Change in fMRI-based rsFC between PMC/SMA seed (asterisk) and dorsal visual stream. * * * p <.001.

Fig. 6.

Correlation across fMRI and EEG measures: Correlations between regional β-ERD distributions and regional fMRI activation patterns. r = 0.81, p <.0001.

During Visual-cortex stimulation, a significant reduction in activity was observed particularly over PMC/SMA, dorsal (Fig. 5A, arrow) and ventral visual cortex. When these regions were used as seeds in a gPPI analysis, a significant negative correlation was observed between both motor and dorsal visual cortex and SMA at baseline. Visual tDCS significantly altered the pattern of dorsal-SMA interaction, as well as local interaction within the dorsal visual region (Fig. 5B). Furthermore, a persistent change in SMA to dorsal visual interaction was observed even following cessation of tDCS (Fig. 5C).

SRT results

SRT behavior

RT was also assessed in a pure SRT condition, in which subjects pressed with a single finger whenever they saw a centrally located stimulus. Data were collected before, during and after stimulation. ERP were obtained only during stimulation. In all cases, data were best fit by a single Gaussian distribution (all goodness of fit R2 > 0.95) (Fig. 7A).

Fig. 7.

Simple Reaction Time (SRT) response: A. Single trial RT distributions: showing single Gaussian distribution under all stimulation conditions. B. Mean RT by tDCS stimulation type: * * p <.01; * * * p <.001 vs sham. C. Correlation: between RT and β-activity at PMC/SMA region under sham and motor anodal stimulation conditions. D. Coherence: significant differential coherence measures when comparing SRT vs fast SRTT trials, SRT vs slow SRTT trials, and random vs SRT. * p <.05.

At baseline, the difference in mean RT across tDCS conditions was not significant (F3,64 = 2.41, p=.08). By contrast, during stimulation, the mean RT varied significantly across the tDCS conditions (F3,64 = 5.98, p=.0012), with significant pairwise reductions for both Visual (F1,32 = 5.81, p=.02) and Motor-anodal (F1,32 = 14.8, p=.0005) stimulation vs. Sham (Fig. 7B).

Similar effects were obtained by ANOVA across all trials and conditions, which showed a significant main effect of condition (F3,14,568 = 23.3, p <.0001), with significant post-hoc differences between both Visual and Motor-anodal stimulation and sham during stimulation. Despite significance, the degree of RT reduction was relatively modest (~15 ms).

SRT neurophysiology

In β-analyses, the main effect of tDCS was not significant (F3,45,204 = 2.20, p=.086) although a significant Region X Condition effect was observed (F6,90,406 = 3.93, p=.004). When analyses were performed by region, tDCS effects were significant only for Motor cortex (F3,45,204 = 6.76, p <.0001), with Motor-cathodal stimulation decreasing β-activity vs. sham (p=.008).

In a mixed-model regression of RT vs. regional β-activity across all stimulation conditions, there were overall main effects of PMC/SMA (F1,45,176 = 4.94, p=.026) β-activity and tDCS stimulation-type (F3,45,176 = 19.7, p <.0001), and significant interactions of tDCS with PMC/SMA (F1,45,176 = 8.25, p <.0001), Visual (F1,45,176 = 5.42, p=.001), PMC/SMA × Motor (F1,45,176 = 8.02, p <.0001) and Motor × Visual (F1,45,176 = 825, p <.0001) sources.

Separate analyses were therefore performed by tDCS condition. Under Sham, single-trial RT correlated primarily with single-trial β-activity in Visual cortex (F1,12,365 = 6.70, p=.01) and the interaction of PMC/SMA and Motor (F1,12,365 = 10.1, p=.001) sources, while the relationship with β-activity in PMC/SMA was not significant (F1,12,365 = 3.46, p=.063). By contrast, during anodal stimulation, a highly significant correlation was observed between RT and single-trial β-activity in PMC/SMA (F1,11,965 = 17.1, p <.0001) (Fig. 7C) and with the interaction between PMC/SMA x Motor sources (F1,11,965 = 12.5, p <.0001). No significant differences were observed between Sham and Motor-anodal stimulation for the lowest 2 hexiles (p >.2), but highly significant differences were observed for each of highest 4 hexiles of β-ERD (all p <.002).

Comparison of coherence between SRT and fast trials of SRTT showed no significant differences across the PMC/SMA - Motor - Visual sources, supporting the hypothesis that fast SRTT trials represent a functional shift to an SRT-like process. Comparison of coherence between SRT vs slow trials of SRTT showed increased coherence across PMC/SMA and Motor (0.18 ± 0.04, p=.03) and across Motor and Visual (0.09 ± 0.03, p=.04) cortices. Finally, Comparison of coherence between random trials of SRTT vs SRT showed increased coherence across PMC/SMA and Motor (0.21 ±.05, p=.04) cortices (Fig. 7D).

Discussion

tDCS is a non-invasive brain stimulation approach in which low level (1–4 mA) currents are applied across the scalp (Bikson et al., 2006). Although the technique is widely used in both basic research and clinical treatment studies, predictability and reproducibility of effect remain major issues (Bikson et al., 2018; Parkin et al., 2019). Here, we evaluated the sensitivity of EEG- and fMRI-based connectivity measures to tDCS effects and motor learning in the SRTT task. Our main findings are threefold.

First, we demonstrate that in the SRTT motor learning is associated with a shift in response mode from slow, reactive CRT-like response to fast, predictive SRT-like responses, similar to what has been observed in other motor learning tasks (Liebrand et al., 2017; Yu et al., 2019). Moreover, we demonstrate that these response types show differential connectivity patterns, such that slow responses are associated with increased β-coherence between dorsal visual cortex and SMA/PMC relative to fast responses. This finding is consistent with a model of behavioral plasticity in the SRTT in which, during fast response trials, the subject has already identified where the stimulus will appear, and only needs to use the visual stimulus to determine timing of the response. In such cases, the response depends primarily on direct interaction between dorsal visual and motor cortices, with limited involvement of PMC/SMA.

By contrast, during slow responses, visual cortex must communicate both location and timing of stimulus information to PMC/SMA, which then determines the response type, and then transmits to motor cortex. The longer response latency in the slow trials thus reflects the differential involvement of PMC/SMA. Consistent with this hypothesis, in the SRTT “random” condition, in which only slow responses are possible, connectivity was significantly higher between SMA and motor cortex than it was during the fixed condition. Consequently, the shift from slow to fast responses may be seen as resulting from competition between the visual → motor and visual → PMC/SMA projections, with the fast responses reflecting a “win” for the more direct pathway.

Second, we demonstrate that effects of tDCS on motor learning in the SRTT are related to a modulation of the shift between fast and slow responses, rather than a shift in mean latency of either response type, and that this shift is related to alterations in functional connectivity patterns as reflected in both EEG β-coherence and fMRI gPPI measures. Prior studies of cortical stimulation in the SRTT have focused primarily on frontal, premotor and motor cortex (Nitsche et al., 2003; Kang and Paik, 2011; Robertson, 2005; Robertson et al., 2001). By contrast, learning of the motor sequence is thought to be encoded within the dorsal visual stream, also termed the “perception for action” system (Keele et al., 2003; O’Regan and Noe, 2001). Stimulation over dorsal visual stream would thus plausibly improve plasticity in the SRTT by increasing the efficiency of motor encoding

Third, we further document that the “anodal excitatory, cathodal inhibitory” (AECI) model of predicting tDCS effect does not hold for motor learning/plasticity, although it may hold for simpler types of motor response tasks. The AECI model is based upon the concept that anodal stimulation produces net depolarization in the apical dendrite of the superficial cortical pyramidal neurons due to its asymmetric orientation. This depolarization then increases cortical excitability, leading to decreases in both TMS motor threshold (i.e., less energy needed to activate a response) and shorter RT in a simple finger tapping task. Cathodal stimulation has the opposite effect (e.g. (Nitsche and Paulus, 2000; Nitsche et al., 2003; Antal, 2004)). However, the degree of change in RT noted as a consequence of anodal stimulation in a simple motor response task (e.g., ~15 ms, Fig. 7), is very small compared to the much larger changes in RT seen in motor learning tasks such as the SRTT (e.g., ~300 ms, Fig. 2), suggesting alternate mechanisms.

An overarching goal of this study was to provide side-by-side evaluation of EEG- and fMRI-based measures that could serve as target engagement biomarkers during clinical treatment studies. To date, tDCS has been used extensively in healthy volunteer subjects in order to assess effects on various cognitive functions. On the other hand, it has been used to only a limited degree in clinical trials for treatment of specific neuropsychiatric conditions, and its utility for these indications has been severely hampered by lack of target engagement biomarkers. Specifically, NIMH-supported early-stage testing of device-based intervention studies require demonstration of a treatment effect on a target engagement biomarker before initiating a clinical trial (Bikson et al., 2018). In the absence of such biomarkers, it remains unknown whether negative clinical trials were unsuccessful because the underlying theory was incorrect, or because the treatment simply did not engage the intended brain circuit.

In our EEG studies, we evaluated effects on ongoing EEG activity, surface-based ERP (Fig. 3), source-space β ERD (Fig. 4A) and source-space connectivity (β-coherence) (Fig. 4B,C,D). No significant effects were observed in ongoing EEG, consistent with observations that such changes only occur with voltage gradients > 1 mV/mm (Vöröslakos et al., 2018). Surface ERP measures were sensitive to tDCS effects especially when analyzed in the frequency domain (Fig. 3D), but the effects were too limited for detailed mechanistic analysis.

More robust effects were detected when ERP responses were mapped into source space using a “beamformer” approach. When this was done, a significant relationship was observed between changes in motor learning and tDCS-induced network-level alterations using either β-coherence (EEG) or gPPI (fMRI). β-coherence measures, in particular, were sensitive to the differential network involvement in fast vs. slow trials and across different task versions (Fig. 4B,C,D). Significant effects were observed for both of the tDCS conditions that enhanced motor learning (motor-cathodal, visual-cathodal), but not for the motor-anodal condition that was ineffective during the SRTT (Fig. 4B, lower right).

Importantly, both motor cathodal and visual cathodal stimulations appeared to function by altering different connectivities within the distributed motor learning network. Thus, whereas cathodal tDCS over motor cortex primarily down-modulated connectivity between dorsal visual stream and PMC/SMA, potentially reflecting increased processing efficiency within that pathway, cathodal tDCS over dorsal visual cortex primarily downregulated β-coherence between dorsal visual and motor cortex, suggesting differential mechanism of action.

In fMRI, we investigated both BOLD activation over the course of the task and effective connectivity as assessed using gPPI. In these studies, significant reductions in BOLD response were observed over both dorsal visual and motor cortex with visual-cathodal stimulation (Fig. 5A). Moreover, across brain regions we observed a highly significant (p <.0001) correlation between ERP and fMRI activation patterns, providing convergent validity between the two approaches (Fig. 6). In gPPI studies, the effects of dorsal visual stimulation were associated with alterations in effective connectivity within dorsal stream itself, as well as within the 3-way SMA, dorsal and dorsal visual network.

Finally, in pre/post rsfMRI, the observed persistent improvement in motor learning was associated with a significant reduction in dorsal visual to motor connectivity (Fig. 5C). Overall, these findings especially support the use of connectivity measures in assessing tDCS target engagement during learning tasks that require distributed neural networks. EEG measures are far simpler to obtain concurrently with tDCS and show strong sensitivity to tDCS effects. Nevertheless, fMRI-based connectivity measures are feasible, provide improved spatial localization and may be used in further refining electrode placement or in guiding high-definition tDCS approaches.

Our findings are consistent with an emerging literature demonstrating network-level effects (Mehrkanoon et al., 2016; Breakspear, 2017) of tDCS (Heth and Lavidor, 2015). For example, Polanía et al. (2018, 2011 a, Polanía et al., 2012, 2011 b) demonstrated significant alterations in EEG and fMRI connectivity during a repetitive motor task, although they did not investigate the interrelationship between connectivity and behavioral changes. Subsequently, it has been demonstrated that effects of tDCS on BOLD activation were dissociable from its local excitatory effects as assessed using TMS, suggesting diverse physiological mechanisms (Antal et al., 2011 a; Lopez-Alonso et al., 2015). Beneficial effects of tDCS on disorders such as tinnitus, depression and dyslexia may also be attributable to network-level effects (To et al., 2018).

Our behavioral results are also somewhat more robust than those observed in other tDCS studies, (rev. in Buch et al. (2017)) even for tDCS applied over motor cortex. This may be attributable to several factors. First, we used a current of 2 mA, as opposed to prior studies that primarily used stimulation currents of 1 mA. One prior study (Cuypers et al., 2013) found numerically, but not significantly, greater change in RT following 1.5 vs. 1.0 mA stimulation and obtained improvements in RT similar in magnitude to those observed here. Use of the higher currents in the present study was enabled by recent safety recommendations that support higher current levels (Nitsche and Bikson, 2017). The low discomfort levels associated with tDCS stimulation in the present study further support use of 2 mA stimulation in future tDCS studies.

A limitation of this study is the relatively small sample size, especially for the imaging components (n = 8). Nevertheless, our sample size is similar to that of other studies that investigated tDCS effects on the SRTT (rev. in Buch et al. (2017)) and sufficient to detect significant tDCS effects. Moreover, we note that this is the first study of which we are aware to combine tDCS treatment with either ERP or fMRI in the context of the SRTT. In addition, order of the Visual stimulation condition was not counterbalanced across subjects, which would potentially affect the across-group comparisons, but would not account for differential patterns of effect associated with the different conditions. Moreover, the Visual stimulation effect remained significant even following control for stimulation order. Nevertheless, future studies using a counterbalanced design are required to replicate the findings.

In addition to its specific findings regarding tDCS modulation of motor learning in the SRTT in healthy individuals, the present study may also have direct implications for use of network-targeted tDCS for treatment of motor manifestations of neuropsychiatric disorders. Thus, it has become increasingly appreciated that “whole brain” diseases such as schizophrenia have a prominent motor component including impairments in TMS-induced motor plasticity (Mehta et al., 2019) that may significantly contribute to poor treatment outcome (Green et al., 2004; Morrens et al., 2014). To this end, a motor domain has recently been added to the Research Domain Criteria (RDoC) matrix of the NIMH, facilitating future research (Bernard and Mittal, 2015). In schizophrenia, reductions in motor dexterity as measured by tests such as the “grooved pegboard” are as robust as those observed in other cognitive domains, and strongly predictive of outcome (Revheim et al., 2006; Lehoux, 2003; Lin et al., 2015; Dickson et al., 2020; Galderisi et al., 2009). Future studies investigating effects of tDCS on motor learning in schizophrenia and other neuropsychiatric disorders thus appear warranted as well.

Conclusions

We provide the first demonstration that 1) motor learning in the SRTT reflects a shift from slow, CRT-like responses to fast, SRT-like responses, 2) the shift in response type reflects a shift in relative connectivity of dorsal visual cortex to PMC/SMA vs. motor cortex, 3) tDCS modulates these alterations, and 4) beneficial effects may be obtained by specific stimulation approaches targeting differential components of the circuit. We further demonstrate that EEG and fMRI approaches provide complementary and convergent information regarding underlying mechanisms. Finally, we demonstrate that source-level β-coherence measures obtained concurrent with tDCS may serve as an effective “target engagement” biomarker for motor learning, potentially allowing greater application of this approach to both neurological and psychiatric disorders associated with prominent impairments in motor plasticity.