Abstract
Conscious intention is a fundamental aspect of the human experience. Despite long-standing interest in the basis and implications of intention, its underlying neurobiological mechanisms remain poorly understood. Using high-definition transcranial DC stimulation (tDCS), we observed that enhancing spontaneous neuronal excitability in both the angular gyrus and the primary motor cortex caused the reported time of conscious movement intention to be ∼60–70 ms earlier. Slow brain waves recorded ∼2–3 s before movement onset, as well as hundreds of milliseconds after movement onset, independently correlated with the modulation of conscious intention by brain stimulation. These brain activities together accounted for 81% of interindividual variability in the modulation of movement intention by brain stimulation. A computational model using coupled leaky integrator units with biophysically plausible assumptions about the effect of tDCS captured the effects of stimulation on both neural activity and behavior. These results reveal a temporally extended brain process underlying conscious movement intention that spans seconds around movement commencement.
Introduction
The two cardinal dimensions of our day-to-day conscious experience are perception and volition (Gray et al., 2007). Volition refers to the sense that one is freely choosing one's actions. It provides the foundation for an individual to attribute agency to the self, and for society to attribute responsibility to an individual. A distorted sense of volition is a hallmark of many neurological and psychiatric illnesses such as alien hand syndrome, psychogenic movement disorders, and schizophrenia (Scepkowski and Cronin-Golomb, 2003; Jeannerod, 2009; Hallett et al., 2012). Thus, understanding the biological basis of volition has paramount scientific, philosophical, legal, and clinical value.
In contrast to the folk psychological view that volition is constituted by conscious intention-causing action, current neuroscientific data favor the view that the sense of volition occurs when conscious intention accompanies self-generated action, and an apparent mental causation is drawn between intention and action (Wegner, 2003). Supporting this view, conscious intention and movement can be dissociated under brain stimulation or pathological states (Fried et al., 1991; Blakemore et al., 2002; Hallett, 2007; Desmurget et al., 2009). In an influential study, Libet et al. (1983) asked subjects to make spontaneous finger movements and report the time at which they first became aware of their intention (W) to move (W-time) using a fast revolving clock. This time was determined to precede movement onset by ∼200 ms, much later (∼1–2 s) than the unconscious initiation of movement indexed by the motor readiness potential [bereitschaftspotential (BP)], supporting the view that conscious intention accompanies but does not cause movement. Despite intense interest in this phenomenon (Lau et al., 2004; Sirigu et al., 2004; Haggard, 2008; Soon et al., 2008; Fried et al., 2011), the neurobiological underpinning of conscious intention remains elusive. One puzzling set of experiments showed that an event happening tens to 200 ms after movement onset could still alter W-times reported to be ∼200 ms before movement (Lau et al., 2007; Banks and Isham, 2009), yet the brain mechanisms mediating such an effect have never been elucidated. Moreover, it remains debated whether conscious movement intention is generated by an instantaneous event (e.g., when the underlying neural activity crosses a threshold; Fried et al., 2011), or by a temporally extended process spanning at least hundreds of milliseconds (Lau et al., 2007; Matsuhashi and Hallett, 2008).
We used high-definition transcranial DC stimulation (HD-tDCS; Datta et al., 2009) to modulate spontaneous neural activity at selected nodes in the motor–premotor–parietal circuit and investigated the resulting effect on conscious movement intention in the context of self-generated movement. It has been suggested that tDCS is an ideal tool for modulating slow brain potentials such as the BP (Elbert et al., 1981). Surface-anodal tDCS has been shown to enhance neural excitability by inducing tonic depolarization of resting membrane potential and increased spontaneous neuronal firing rate (Nitsche et al., 2008). The aftereffect of HD-tDCS, mediated by BDNF-dependent long-lasting synaptic potentiation (Fritsch et al., 2010), has been shown to last for >2 h (Kuo et al., 2013).
Materials and Methods
Subjects
The study was approved by the Institutional Review Board of the National Institute of Neurological Disorders and Stroke. Twelve right-handed healthy volunteers (seven female) with normal or corrected-to-normal vision participated in the study. The study used a within-subject design, with each subject undergoing six different tDCS treatments (see below). Their mean age was 25.6 years (age range, 23–38 years). All subjects provided written informed consent.
Overall experimental design
Each subject participated in six experimental sessions. In each session, subjects received anodal or sham tDCS, and immediately afterward performed the behavioral task under electroencephalography (EEG) recording. The behavioral EEG portion of the experiment was conducted within a window of ∼1.5 h after the end of tDCS stimulation. Due to the time required to prepare EEG electrodes, the start of the experiment ranged from 18 to 33 min (mean, 23.6 min) after the end of tDCS stimulation, and the end of the experiment ranged from 60 to 100 min (mean, 70.6 min) after the end of tDCS stimulation. Three different brain regions were targeted, each of which was stimulated with anodal and sham tDCS in different sessions. Sessions were scheduled on separate days, with adjacent sessions spaced at least 7 d apart. Previous studies (Kuo et al., 2013) have found that the effect of a single session of HD-tDCS (at 2 mA for 10 min) lasts for >2 h after stimulation ends. Anodal and sham stimulation to the same brain region was always performed in adjacent sessions to mitigate potential long-term drift in subjects' performance over multiple sessions, while the sequence of stimulated regions and the order between anodal and sham stimulation were counterbalanced across subjects. In all analyses, sham data from the three sessions were pooled together, unless otherwise noted.
Neurotargeting
In all subjects, anatomical magnetic resonance imaging (MRI) data were acquired on a General Electric 3 T scanner with an eight-channel head coil, using an MP-RAGE sequence with a resolution of 1 × 1 × 1 mm3.
The following three brain regions were selected for stimulation: supplementary motor area (SMA)/pre-SMA, contralateral (left) primary motor cortex (M1), and contralateral (left) angular gyrus (AG). The locations of these brain regions were determined based on coordinates published in prior studies (Table 1, Fig. 1B). For SMA/pre-SMA, coordinates were determined as the average between SMA and pre-SMA to stimulate both regions equally, since the spatial focality of tDCS precludes stimulation of each region separately.
Anatomical locations for stimulated brain regions
Stimulation sites were determined based on each individual subject's anatomical MRI, using the Brainsight neuronavigation system (Rogue Research Inc.). Each subject's MRI data were registered to the standard MNI space using custom AFNI scripts for nonlinear transform (Cox, 1996). After the targeted coordinates were applied to MNI space, the MRI was transformed back into the individual subject's space for targeting. A mark was made on the skin at the scalp location nearest to the target region for placement of the tDCS anode.
tDCS
tDCS was administered under the supervision of a credentialed physician or nurse practitioner. Current was generated by a neuroConn DC Stimulator Plus channeled through a 4 × 1 Multichannel Stimulation Device (Soterix Medical). A 4 × 1 montage consisting of five sintered Ag/AgCl ring electrodes (Stens; Minhas et al., 2010) was used, consisting of one anode directly over the stimulation site surrounded by four equally spaced return electrodes (i.e., cathodes) at a radius of 5 cm from the anode (Fig. 1A). The electrodes were held in place in plastic electrode holders in a fitted cap (EASYCAP). The electrode holders were filled with SignaGel, creating a gel contact of ∼4 cm2 per electrode (peak current density at the anode, ∼0.5 mA/cm2; peak current density at the cathodes, ∼0.125 mA/cm2).
A 4 × 1 ring HD-tDCS has been used to replicate findings of conventional tDCS, which uses saline-soaked sponge electrodes with a bipolar montage. HD-tDCS has been shown to have better spatial focality, larger effect on cortical excitability, and longer aftereffect than the conventional bipolar tDCS (Datta et al., 2009; Caparelli-Daquer et al., 2012; Kuo et al., 2013). In addition, HD-tDCS has been shown to be safe and tolerable at 2.0 mA for 20 min, and allows for effective sham (Borckardt et al., 2012; Villamar et al., 2013).
Anodal stimulation consisted of a linear ramp-up period of 30 s at the beginning, followed by 20 min of sustained stimulation at 2.0 mA, and, last, a linear ramp-down period of 30 s. Sham stimulation consisted of a 30 s ramp-up period to induce a tingling sensation, after which the current was turned off abruptly, and the subject waited for 20.5 min. During both anodal and sham tDCS, subjects were asked to relax and remain awake. For the effectiveness of sham stimulation, see Sham effectiveness in Results.
EEG recording
Subjects were seated in a dimly lit, electromagnetic interference-shielded, soundproof room. Stimuli were presented on an LCD monitor (resolution, 1920 × 1080; refresh rate, 120 Hz). Sixty-four Ag/AgCl actiCAP EEG electrodes (Brain Products GmbH) were placed according to the international 10–20 system. Sixteen additional electrodes were placed over central frontal areas to allow better recording of the motor readiness potential (Fig. 2). Two reference electrodes were placed on the left and right mastoid. The following four electro-oculogram (EOG) electrodes were placed: two above and below the right eye; and two at the outer corner of both eyes. Two electromyogram (EMG) electrodes were placed 2 cm apart over the right index finger flexor muscle (flexor digitorum superficialis). All electrodes were prepared with 10% chloride ABRALYT HiCl abrasive electrolyte-gel (EASYCAP). EEG, EOG, and EMG data were collected using the BrainAmpDC system (Brain Products GmbH) in DC recording mode with a sampling rate of 500 Hz.
Behavioral task
Self-paced movement task.
Subjects performed a self-paced movement task in which they were asked to click a mouse with their right index finger at a time of their choosing (Fig. 1C). Each trial started with the appearance of a clock in the center of the screen, with a black dot positioned at a random location around its perimeter. After a 1 s delay, the black dot changed to red and began to rotate around the clock with a period of 2.56 s per rotation. Subjects were asked to allow the dot to make at least one full rotation before clicking but to allow the urge to move appear spontaneously on its own, without any preplanning of when to move, and then to click the mouse immediately upon feeling the urge. After clicking, the dot continued to rotate for a random period between one to two full rotations. Subjects were then prompted to report the timing of movement (M-time) or awareness of W by using the mouse to precisely position the red dot at the same position as the onset of movement or awareness of intention. After a fixed intertrial interval (ITI) of 4 s, the next trial began. In M-trials, subjects were asked to report the location of the red dot on the clock when they began to move their finger. In W-trials, subjects were asked to report the location of the red dot on the clock when they first became aware of their intention to move.
Control task.
To control for the precision in using the clock, subjects also performed an auditory timing task [sound (S) trials; Fig. 1D]. This task approximately matched the format of the self-paced movement task except that rather than making a spontaneous movement, subjects listened for a brief auditory tone (1 kHz) presented for 100 ms. The interval between the onset of the red rotating dot, and the onset of the tone was randomly chosen between 2.56 and 10.24 s (i.e., one and four full rotations of the dot). After the tone ended, the dot continued to rotate for a random period between one and two full rotations. Subjects were then asked to report the time on the clock when they heard the tone.
In the original study by Libet et al. (1983), the reported time of an S-stimulus delivered at random times was used to “correct” for potential error in the timing report during M- and W-trials by subtracting the reported S-time from M- and W-times on a session-by-session basis. However, although the S-trials share with M- and W-trials the process of monitoring the revolving clock and reporting the timing of an event based on it, monitoring an unpredictable external event (skin or auditory stimulus) is unlikely to be a subset of processes involved in monitoring an internal event (movement intention or motor awareness). Thus, in line with recent prior studies (Sirigu et al., 2004; Lau et al., 2007), we opted to use the S-trials as a separate control condition instead of correcting for M- or W-times based on S-times.
Task structure.
Subjects first completed five practice trials of each task type, followed by two blocks of M-trials, two blocks of W-trials, and one block of S-trials. Each block included 30 trials. The block order was staggered and counterbalanced across subjects. The task paradigm was implemented in Python and presented using the OpenSesame software package (Mathôt et al., 2012).
EEG signal preprocessing
EEG data preprocessing was performed using BrainVision Analyzer 2 software (Brain Products GmbH), Matlab (MathWorks), and EEGLAB (Brunner et al., 2013).
First, scalp and eye electrodes were referenced to digitally linked mastoid reference electrodes. Continuous EEG and EOG data were filtered with a bandpass filter between 0.01 and 150 Hz (24 dB/oct) and a 60 Hz notch filter. All filters were applied off-line using a symmetrical Butterworth filter with zero phase shift. Second, blinking and eye movement artifacts were removed using an independent component analysis (ICA) applied to all EEG and EOG electrodes. ICA was computed using the InfoMax Extended algorithm implemented in Analyzer 2 with a maximum of 512 steps and a convergence bound of 1 × 10−7. Next, EEG signals were epoched from 4 s before to 2 s after the button press in M- and W-trials, and from 4 s before to 2 s after the tone onset in S-trials. Baseline correction was performed using the first 500 ms of each segment (i.e., 4000–3500 ms before button press), and linear detrending was applied to remove slow drift. Last, artifact rejection was performed by removing all trials with ≥1000 μV fluctuations or activity >5 SDs from the mean using EEGLAB.
EMG electrodes were referenced in a bipolar montage. A 0.01–50 Hz bandpass filter (24 dB/oct) and a 60 Hz notch filter were applied to the EMG signal, which was then epoched similarly to the EEG data. Movement onsets were identified using an algorithm that detected the first data point with both a positive second derivative and an activity level >4 SDs above or below the 500 ms baseline period at the beginning of each segment. The automatic algorithm was validated by visual inspection of every trial. Averaged across all tDCS conditions, EMG onset was at 158 ± 3 ms (mean ± SEM across subjects) before the onset of button press in W-trials, and 157 ± 3 ms before button press onset in M-trials. None of the anodal stimulation conditions significantly altered the EMG onset compared with sham (all p > 0.9, random-effects population t test).
Event-related potential analysis
For event-related potential (ERP) analysis, an additional <15 Hz low-pass filter was applied to epoched data using a Hamming windowed sinc FIR filter in EEGLAB.
Motor readiness potential.
A pool of 17 electrodes centered around C1 was selected based on the strength of the BP signal in the sham condition in a 300-ms-long window right before button press. These electrodes were Cz, C1, C3, C5, CCP1h, CCP3h, CCP5h, CP1, CP3, FCC1h, FCC3h, FCC5h, FCz, FC1, FC3, FFC1h, and FFC3h (Fig. 2). Electrode nomenclature follows the 128-channel high-density EEG based on the extended International 10–20 system (Oostenveld and Praamstra, 2001).
Lateralized readiness potential.
Traditionally, lateralized readiness potential (LRP) is defined as a double subtraction, as follows: C3–C4 difference potential for trials with right-hand responses (RH) is subtracted from C3–C4 difference potential for trials with left-hand (LH) responses [LRP = (C3–C4)RH−(C3–C4)LH] (Jahanshahi and Hallett, 2003). We calculated the left component of the LRP under right-hand responses [(C3–C4)RH]. Previous work has shown that this is a fair measure of lateralization and that, due to deviations in brain symmetry, the double subtraction method may introduce error and does not provide significant benefit over a single subtraction (Oostenveld et al., 2003).
Parietal cortical activity.
Using a Brainsight neurotargeting device on subjects wearing a cap with standard 10–20 locations, we determined that a pool of electrodes including P1, P3, and PO3 were directly over the left angular gyrus (as defined by MNI coordinates in Table 1; Fig. 2). This selection of electrodes is consistent with previous literature, which often used the P3 electrode to index activity of the inferior parietal lobule (Herwig et al., 2003).
Statistical tests
For behavioral data, M- and W-times were referenced to the onset of button press, and S-times were referenced to the onset of auditory tone. A two-sample random-effects t test including both within-subject across-trial variances and across-subject variances (Borenstein et al., 2007) was used to compare behavioral results between anodal and sham condition. Two different analyses were conducted. In the first, behavioral data under anodal stimulation of each brain region were compared with sham data pooled across all three regions. In the second, behavioral data under anodal stimulation of each brain region were compared with sham stimulation of the same brain region.
Additionally, repeated-measures ANOVA with factors of stimulation site (AG, M1, SMA) and treatment (anodal, sham) was conducted to provide an overall assessment of how the experimental manipulations affected behavioral measures. This ANOVA design was applied to the mean and SD of W-time and M-time, as well as to the mean and SD of the time from trial onset to button press (“click time”). The purpose of these analyses was to investigate the specificity of the behavioral effect of tDCS.
For ERP analyses, EEG signal amplitude was compared between anodal (of each brain region) and sham (pooled across three regions) conditions using a two-sample t test without assuming equal variances on every time point in the segmented epoch. Results were corrected for multiple comparisons using Bonferroni correction with a correction factor derived from Bartlett's theory (Jenkins and Watts, 1998; Vincent et al., 2007; Van Dijk et al., 2010; He and Zempel, 2013). This correction factor takes into account the amount of autocorrelation in the data, which was used to derive the true degree of freedom and the number of independent tests.
To assess the relation between tDCS influences on W-time and brain activity, we computed the across-subject Pearson correlation between the effect of tDCS on W-time and its effect on ERP amplitude. ERP was averaged across trials for each subject. The effect of tDCS on both behavior and brain activity was calculated as the difference between the anodal and sham conditions. A search algorithm was used to identify time windows in which the effect of tDCS on ERP correlated significantly with its effect on W-time, which used 300-ms-long windows that advanced in 50 ms steps. The results were corrected for multiple comparisons using Bonferroni correction with the Bartlett correction factor derived for ERP data smoothed with 300-ms-long, 50 ms-step moving averaging windows. Final windows selected for windowed tests were determined by the first and last time point in groups of adjacent significant search windows. To exclude potential outliers, windowed tests were further confirmed by robust linear regression analyses.
Computational model
To capture the main features of the EEG and behavioral data as well as their modulation by tDCS, we constructed a simple yet biophysically plausible model following the spirit of the mean-field approach (Wong and Wang, 2006). The model implements three interacting leaky integrator units (beim Graben et al., 2007), corresponding to SMA/pre-SMA, M1/premotor cortex (PMC), and AG (see Fig. 6B). We take activity from the M1/PMC node as a proxy for the LRP, activity from the AG node as a proxy for parietal activity, and activity from the SMA/pre-SMA node as a proxy for the BP. The model makes the following simplifications: (1) it explicitly models only SMA/pre-SMA, M1/PMC, and AG activity, thus, it does not explicitly model other potentially relevant regions, such as the basal ganglia, primary sensory cortex, and the ipsilateral M1/PMC; (2) given the limited spatial resolution of tDCS and EEG, the model does not differentiate between SMA and pre-SMA, or between M1 and PMC; (3) it posits a simplified pattern of connectivity among these regions; and (4) in definitions of specific parameter values, the general features of the parameters such as their sign and time course were determined a priori by theoretical and empirical considerations, while the specific numerical values are set to capture relevant features of the empirical data without an exhaustive parameter search.
In general, the rate of change of a leaky integrator x as a function of time is described as follows:
where τ is a time constant controlling how rapidly the activity decays over time, and y(t) is the input to x at time t. For simulation purposes, we used a discrete time axis incrementing in steps of Δt = 1 ms. Let xSMA(t), xM1(t), and xAG(t) denote the activity at time t for the leaky integrator units corresponding to SMA/pre-SMA, M1/PMC, and AG, respectively. For each unit, the time course of activity was evaluated over the interval t = [−10 s, 2 s], where t = 0 represents the time of button press. The activity at the initial time point t = −10 s for all units was set to 0. From then on, activity was calculated iteratively at each time step. For all units, we set τ = 1000 ms to produce simulated neural activity with rates of accumulation before button press, and rates of decay after button press, comparable to those observed in the empirical ERP. Mechanistically, such a long time constant could be implemented by recurrent neural networks within each brain region (Wang, 1999; Major and Tank, 2004).
Although our empirical ERP results were assessed only within a time window of [−4 s, 2 s] around movement, we simulated a much longer time window in the model. This long time window was used to capture the tDCS modulation of baseline activity, as specified by the parameters bM1 stim and bAG stim described below. During t = [−10 s, −4 s], these baseline modulation parameters integrate to a stable asymptote, corresponding to a shift in baseline activity. To compare simulated neural activity with empirical ERPs, which were baseline corrected (see EEG signal preprocessing), simulated activity in the time window [−4 s, 2 s] was also baseline corrected using the activity at −4 s.
The input term y(t) differs for each unit, depending on the connectivity of that unit to the other units and external inputs (see Fig. 6B). In the simulation of the sham condition, y(t) was defined for each unit as follows.
Input to SMA node.
The input function for SMA is as follows:
where ISMA(t) defines the external input received by SMA, and followed:
to produce slow ramping up of activity in SMA in the 4 s period preceding button press, corresponding to the empirically measured readiness potential. Input turns off after t = 0, allowing SMA activity to decay to baseline following button press. We note that the two-step input function assumed for SMA and AG (see below) is one potential mechanism for generating activity time courses that resemble empirical ERPs. The shape of the ERP under the sham condition could also be produced by a leaky integrator model coupled with a threshold mechanism (Schurger et al., 2012), or nonlinear dynamics within each brain region. The purpose of the model is not to propose detailed mechanistic account for the generation of ERPs under the sham condition; rather, our main focus is on explaining tDCS modulation of ERPs and W-time, as well as the correlation between them, by using biologically based assumptions about the effects of tDCS.
Input to AG node.
The input function for AG is defined as follows:
where IAG(t) defines the external input. AG also receives input from M1, weighted by a factor wM1→AG(t).
IAG(t) followed:
to reproduce the slow ramping-up activity in the parietal ERP, followed by a return to baseline after button press at t = 0.
We supposed that input from M1 to AG primarily reflects computations related to forward modeling, in which a copy of the motor command is relayed to parietal areas to assist in the prediction and on-line control of movement (Wolpert et al., 1995; Desmurget et al., 1999; Sirigu et al., 1999; Ogawa et al., 2007). We therefore modeled the input from M1 to AG as occurring in a 500 ms window around button press, as follows:
A direct connection between AG and M1 in the human brain is supported by prior probabilistic tractography results (Caspers et al., 2011). Following previous literature (Wolpert et al., 1995; Wolpert and Flanagan, 2001; Ogawa et al., 2007), we model the input from M1 to AG as an inhibitory influence for the following considerations: the parietal cortex has been suggested to compare the state estimates computed from efference copy with sensory feedback for on-line error correction (Desmurget et al., 1999; Sirigu et al., 1999; Ogawa et al., 2007); a simple and efficient implementation of such a comparison is a subtraction between the two streams, implemented by simultaneous inputs of excitatory sensory feedback and inhibitory efference copy. This framework also accords well with results showing increased activity in AG when the sense of agency is violated, which is presumably due to incomplete cancellation between sensory feedback and efference copy (Farrer et al., 2003; Nahab et al., 2011).
Input to M1 node.
We defined the input function for M1 as follows:
Thus, in this simple model we account for M1 activity solely as a function of inputs from AG and SMA.
The input from both AG and SMA to M1 were set to be excitatory, with the input from AG to M1 extending a longer period from 4 s before to 500 ms after button press, while the input from SMA occurred in a 1-s-long time window immediately before movement (see Fig. 6C). Specifically, we set the following:
These model choices are based on a framework (Desmurget and Sirigu, 2009, 2012) supported by lesion and cortical stimulation observations, which suggests that the inferior parietal lobule is involved in the initial formation of movement intention (Sirigu et al., 2004; Desmurget et al., 2009) and early movement planning (Wheaton et al., 2005), as well as continued forward modeling and movement monitoring (Desmurget et al., 1999; Sirigu et al., 1999), whereas SMA/pre-SMA emits a “go signal” close to movement onset by releasing the inhibitory control exerted on M1 (Fried et al., 1991; Ball et al., 1999). An excitatory influence from AG to M1 in the left hemisphere was also observed in a transcranial magnetic stimulation (TMS) study (Karabanov et al., 2013).
Modeling the effects of tDCS on network dynamics.
Anodal tDCS enhances spontaneous neural excitability and firing rate by inducing tonic depolarization of resting membrane potential (Creutzfeldt et al., 1962; Purpura and McMurtry, 1965); it is also associated with locally reduced GABA levels (Stagg et al., 2009). In addition, DCS paired with low-frequency synaptic activation causes BDNF-dependent long-lasting synaptic potentiation (Fritsch et al., 2010). Accordingly, our model incorporated two mechanisms whereby tDCS altered network dynamics. First, stimulated nodes exhibit an elevated baseline level of activity, reflecting enhanced excitability. Second, the connectivity strengths among nodes are modulated by tDCS.
Formally, the effect of AG anodal stimulation is captured by changes to the input received by AG, as defined by the following:
There are two modifications from the input function for AG under sham stimulation (Eq. 3). First, the additional term bAG stim captures overall enhanced excitability of AG by elevating its baseline activity at t = −4000 ms. We set bAG stim = 2. Second, the input from M1 onto AG is modulated by a constant kM1→AG|AG stim, set to a value of 3, which reflects enhanced synaptic strength induced by tDCS.
Similarly, the effect of M1 anodal stimulation was defined by the following input to M1:
where bM1 stim = 1, kAG→M1|M1 stim = 1.2, and kSMA→M1|M1 stim = 1.1.
Finally, the effect of SMA anodal stimulation was defined by the following input to SMA:
where bSMA stim = 0.5.
Simulating reported W-times under different tDCS conditions.
We next investigated whether the above model could account for the effect of tDCS on W-time. Our empirical ERP analyses suggested that the effects of tDCS on LRP as early as −2850 ms, and on parietal activity as late as 950 ms, correlated with tDCS modulation of W-time (see Fig. 4). Accordingly, we supposed that the genesis of W-time is implemented by a temporally extended process ranging from 2850 ms before to 950 ms after movement, during which a weighted average of visual input from the rotating clock was computed to construct a final estimate of W-time. We characterized the weights used for this temporal averaging as a function of M1 activity to capture the close correlation between LRP and W-time in the empirical findings (see Figs. 4, 5). We further supposed that not all values of M1 activity in the specified time window were weighted evenly, but rather that activity hundreds of milliseconds before button press would likely be weighted more heavily than activity seconds before or hundreds of milliseconds after button press. Finally, to take into account the ∼80 ms lag between visual stimulus and subjective visual perception (Eagleman and Sejnowski, 2000; Dehaene and Changeux, 2011; see Fig. 6A), we supposed that the weight derived from M1 activity for a given time t in fact multiplied a value of (t − 80).
Formally, we defined W-time as follows:
where W-time is a weighted average of (t − 80) ms in the interval t = [−2850 ms, 950 ms]. The weights w(t) were in turn defined as a weighted function of M1 activity over this time window, and normalized such that w(t) summed to 1:
The β(t) weights on M1 activity were defined as a normal probability density function, as follows:
where μ = −330 ms and σ = 800 ms, to capture the intuition that the weights should be highest hundreds of milliseconds before movement (see Fig. 8A). β(t) was fixed across tDCS conditions; thus, the only source of variation in w(t), and hence W-time, came from variation in xM1(t) across conditions. The time courses of xM1(t) and w(t) under sham, AG and M1 stimulation conditions are shown in Figure 8, B and C.
Simulating interindividual variability in the effects of tDCS.
To assess whether the model could capture the empirical correlation between tDCS modulation of W-time and its effects on ERP, we conducted a Monte Carlo simulation. We simulated 1000 subjects. For each subject, the parameters of the model were defined as above, except that Gaussian noise was added to the parameters modulated by tDCS (i.e., b and k terms). For simplicity, the noise terms added to different parameters modulated by the same tDCS condition were fully correlated across subjects to capture the variation in the effective strength of tDCS across subjects. For each b parameter, the SD of the noise was set equal to the value used for that parameter in the main simulation, with the constraint that b > 0. For each k parameter, the SD of the noise was set as the distance between the main simulation value and the lower bound 1, with the constraint that k > 1.
Testing model predictions
The increased spontaneous neural excitability posited by the model (b terms in Eqs. 5 and 6) predicts baseline shifts in the ERP, which are rendered invisible by the standard baseline correction procedure in EEG signal preprocessing (see Fig. 7). To test this model prediction, we additionally computed LRP and parietal ERP without baseline correction under different tDCS conditions (see Fig. 9A–D). It is important to note that this analysis is expected to be noisier, as the results are affected by slow artifacts in DC-EEG recordings such as electrode drifts over time (Rockstroh et al., 1989).
To test the model prediction regarding increased synaptic efficacy (k terms in Eqs. 5 and 6), we computed functional connectivity (FC) between spontaneous EEG signals using the ITI period. Spontaneous BP, LRP, and parietal activity were extracted from the same electrode groups as for ERP analyses (Fig. 2), and a mean signal value was extracted for each 4-s-long ITI. Functional connectivity was computed as a Pearson correlation between BP and LRP, and between LRP and parietal activity, across ITIs. Correlation values were Fisher z-transformed before averaging across subjects (see Fig. 9E,F). In a second FC analysis, spontaneous EEG signals during the ITI period were smoothed with a 200-ms-long, half-overlapping moving-average window, and concatenated across ITIs; and Pearson correlations were computed on concatenated signals. The results showed an identical pattern to those from the first analysis (data not shown).
Results
Twelve healthy young volunteers participated in the experiment, each performing the experiment by Libet et al. (1983) in six different sessions under different conditions of tDCS. In separate blocks of trials, participants judged the time of when they first became aware of their W-time or when they first began their physical movement (i.e., M-time; Fig. 1C). To control for precision in using the clock, subjects also judged the time of a brief auditory tone (i.e., S-time; Fig. 1D). Subjects performed this task immediately after HD-tDCS (Fig. 1A), within a window of <1.5 h after the termination of stimulation current. While performing the task, subjects' brain activity was monitored by high-density EEG recording. The six different tDCS conditions included anodal and sham stimulation of three brain regions: the SMA/pre-SMA, the contralateral M1, and the contralateral AG (Fig. 1B, Table 1). These brain regions were selected based on prior studies of movement initiation, movement intention, and motor awareness (Fried et al., 1991, 2011; Lau et al., 2004, 2007; Sirigu et al., 2004; Desmurget et al., 2009). The stimulation sites were targeted based on each subject's anatomical MRI, using a neuronavigation system (Fig. 1B).
HD-tDCS setup and task paradigm. A, An example of the tDCS ring montage with the central anode over SMA/pre-SMA and return electrodes (cathodes) spaced at a 5 cm radius. B, Stimulation targets marked on the reconstructed scalp and brain surface of one subject. C, Each trial started with the appearance of a clock in the center of the screen, with a dot positioned at a random location around its perimeter, which turned red after 1 s and began to rotate around the clock with a period of 2.56 s per rotation. After at least one full rotation, subjects clicked the mouse with their right hand when they felt a spontaneous urge to do so. After clicking, the dot continued to rotate for a random period between one and two full rotations. The subject was then prompted to report M or W by using the mouse to precisely position the dot. D, Control task in which, instead of making a spontaneous movement, subjects were asked to report the time when they heard a brief auditory tone that lasted for 100 ms and occurred at a random time between 2.56 and 10.24 s after the clock started rotating.
EEG recording and analyses montage. Eighty scalp electrodes were recorded; their locations are shown on the map. Background color shows the current source density map in a 300-ms-long window right before button press, averaged across all sham sessions. The amplitude of the current source density map was used to select electrodes pooled for analyzing BP (shown in red). The LRP was calculated as a subtraction between electrode C3 (contralateral to hand movement) and C4 (ipsilateral to hand movement). Parietal activity was obtained by pooling across three electrodes overlying the AG (shown in yellow).
Sham effectiveness
At the end of each session, subjects were asked to guess whether they received sham or anodal tDCS in that session, and to rate their confidence level in this response on a scale from 1 to 10 (where 1 indicates “I have no idea” and 10 indicates “I am certain”). Subjects were not able to distinguish anodal from sham tDCS significantly above chance level (p = 0.099, proportion test) and reported overall low confidence in their answers (mean ± SEM, 3.83 ± 0.34). Moreover, their confidence report did not correlate with their objective performance. Their reported confidence was slightly (but not significantly) higher when their guess was wrong (4.24 ± 0.53 in incorrect guesses, 3.56 ± 0.45 for correct guesses; p = 0.29, Wilcoxon rank-sum test). Last, subjects' ability to distinguish sham from anodal stimulation did not increase over time, as assessed by Spearman rank correlation between the percentage of accuracy (pooled across subjects) and session number (from 1 to 6; r = 0.37, p = 0.47). These results indicate that subjects did not have explicit or implicit insight into whether they received anodal or sham stimulation in a given session.
Behavioral results
Under sham condition, the reported W-time was on average 188 ± 22 ms (mean ± SEM across subjects) before the recorded button press, and the reported M-time was 73 ± 15 ms before button press (Fig. 3A). Both values are consistent with previous literature (Jahanshahi and Hallett, 2003). Subjects judged the S-time to occur 54 ± 13 ms after the physical onset of the tone. Both AG and M1 anodal stimulation shifted the reported W-time significantly earlier, by 61 ± 27 ms (p = 0.016, random-effects population t test) and 70 ± 35 ms (p < 0.009), respectively. By contrast, SMA/pre-SMA stimulation did not significantly alter W-time. None of the anodal stimulation conditions significantly changed the reported M-time or S-time (Fig. 3A). These results remained the same when each anodal condition was compared with sham stimulation of the same brain region (Fig. 3B). A two-way repeated-measures ANOVA using W-times as the dependent variable found a significant effect of treatment (anodal vs sham, p < 0.05) and a significant interaction effect (treatment × stimulation site, p < 0.03), while the effect of stimulation site was not significant (p = 0.09). Post hoc analyses suggested that the interaction effect was driven by AG and M1 anodal stimulation, both of which significantly shifted W-time earlier compared with the corresponding sham stimulation (AG: p = 0.003; M1: p = 0.004; random-effects population t tests). The same ANOVA design did not reveal any significant effect of stimulation site, treatment, or their interaction on the SD of W-time or M-time (all p > 0.17).
Behavior results. A, Mean reported S-, M-, and W-times under different stimulation conditions. Sham data were pooled across three sessions in each subject. M- and W-times were referenced to the onset of button press; S-times were referenced to the onset of the tone. Error bars represent the SEM across subjects. Reported times were compared between anodal stimulation of each brain region and pooled sham data, using population-level random-effects t tests. B, Same as A, except that sham results are shown separately for different brain regions, and each anodal stimulation condition was compared with sham stimulation of the same brain region. Asterisks indicate significant comparisons.
We additionally investigated whether tDCS modulated the frequency or variability of subjects' button presses. To this end, we determined the time from trial onset to button press in every trial (click time), and computed the mean and SD of click times across all W-trials or all M-trials, under each tDCS condition. The results are subjected to a two-way repeated-measures ANOVA across subjects [factors: stimulation site (SMA, M1, and AG) and treatment (anodal vs sham)]. All of the main effects and interaction effects were not significant (all p > 0.2), suggesting that tDCS did not alter movement statistics.
Early LRP and late parietal activity independently correlate with W-time modulation
What brain mechanisms might contribute to the effects of AG and M1 stimulation on the reported time-of-movement intention? To address this question, we investigated the following three ERPs chosen a priori for their ability to index neural activity in brain regions modulated by tDCS in this study: the BP (Kornhuber and Deecke, 1965; Libet et al., 1983), which has a generator in the SMA and pre-SMA (Ikeda and Shibasaki, 2003); the LRP (Haggard and Eimer, 1999), which has a generator in M1 (Jahanshahi and Hallett, 2003); and parietal activity from a group of electrodes overlying the AG (Fig. 2; see Materials and Methods).
AG stimulation significantly enhanced the LRP ∼1–3 s before button press (p < 0.05, Bonferroni corrected; Fig. 4) without significantly altering the amplitude of BP or parietal activity. Using a search algorithm to identify time windows during which tDCS modulation of LRP significantly correlated with its modulation of W-time (p < 0.05, Bonferroni corrected), we found an early window in the LRP (2.5–2.85 s before button press)—just before the profound elevation of LRP amplitude—correlating significantly with the behavioral modulation (Fig. 4A,B). The larger the effect of AG stimulation on LRP amplitude in this time period, the earlier the subject reported their movement intention relative to sham (Pearson's correlation: r = 0.71, p < 0.01; robust linear regression: p < 0.016). The early LRP activity alone accounted for 50% of across-subject variance (i.e., interindividual variability) in the modulation of W-time by AG stimulation. In addition, the best-fit regression line passed through the origin (Fig. 4B), suggesting that when early LRP amplitude was unchanged from sham condition, the reported W-time was unchanged. Because AG stimulation was over left parietal cortex, whereas the LRP is a measure of activity lateralization over primary motor cortices, the above finding raises the question of whether a similar pattern would hold for activity lateralization over parietal cortex. A control analysis using parietal activity in the same time window suggested that this was not the case (r = −0.11, p > 0.7).
EEG correlates of the effect of AG stimulation on W-time. A, LRP in W-trials under AG anodal and sham stimulation, averaged across subjects. Negativity is plotted upward per EEG convention. Red dots mark the time points where the two curves significantly depart (p < 0.05, Bonferroni corrected). Purple dots mark the center of 300-ms-long windows during which the effect of AG stimulation on LRP amplitude significantly correlated across subjects with its modulation of W-time [p < 0.05, Bonferroni (Bonf.) corrected]. The yellow window indicates the time period selected for the windowed test in B, as determined by the first and last time points in adjacent significant search windows. B, The effect of AG stimulation on LRP amplitude in a window of 2.5–2.85 s before movement (yellow window in A) is plotted against the effect of AG stimulation on W-time for each subject. Both effects are calculated as the difference between anodal and sham conditions. The ordinary least-squares (OLS) linear regression fit is shown as the solid red line. This is supplemented with a robust regression analysis to assess the potential impact of outliers on the correlation results. The robust regression fit to the data (red dashed line) closely matches the OLS fit, suggesting that outliers do not drive the correlation result. Each data point is labeled with the weight assigned to it by the robust regression algorithm. Lower weights indicate data points that contribute less to the robust regression fit. C, D, Same format as in A and B, for parietal activity. The yellow window in C indicates the time period (0–550 ms after movement) selected for the windowed test shown in D. E, F, E is reproduced from C, except that the yellow window indicates the second temporal cluster of significant parietal activity correlations with W-time, in the period of 450–950 ms after movement. The corresponding windowed test is shown in F.
Although the effect of AG stimulation on parietal activity did not reach significance after Bonferroni correction, parietal activity tended to return to baseline faster after movement under AG stimulation (reaching p = 0.02, uncorrected; Fig. 4C). A search algorithm revealed two contiguous clusters of 300-ms-long windows in which the modulation of parietal activity by AG stimulation significantly correlated with its modulation of W-time (p < 0.05, Bonferroni corrected), with one cluster spanning 0–550 ms after movement, and the other spanning 450–950 ms after movement (Fig. 4C,E). In both time windows, the faster parietal activity returned to baseline, the earlier was the reported W-time (0–550 ms: Pearson's r = −0.78, p < 0.003; robust regression, p = 0.003, Fig. 4D; 450–950 ms: Pearson's r = −0.75, p < 0.005; robust regression, p < 0.006; Fig. 4F). Not surprisingly, these two time windows contained highly redundant information. The effect of AG stimulation on parietal activity in the two windows are highly correlated across trials (r = 0.85, p < 0.0005). Moreover, while the first time window alone accounted for 61% of across-subject variance in the modulation of W-time by AG stimulation, a general linear model incorporating both windows had only marginal improvement, explaining 64% of behavioral variance. Hence, in the following analysis, we focus on parietal activity in the first time window (0–550 ms).
Since both early LRP and late parietal activity correlated with the effect of AG stimulation on W-time, an important question is whether they contributed to the behavioral modulation independently. Two pieces of evidence suggest that this is the case. First, the effects of AG stimulation on early (2.5–2.85 s before movement) LRP and late (0–550 ms after movement) parietal activity were not correlated across subjects (p > 0.2). Second, using partial correlation, we found that controlling for the effect of late parietal activity, the early LRP window still significantly correlated with the modulation of W-time (r = 0.71, p < 0.014); similarly, controlling for the effect of early LRP, the late parietal activity still significantly correlated with W-time modulation (r = −0.79, p < 0.005). A general linear model incorporating both early LRP and late parietal activity as predictors suggested that they together accounted for 80.9% of interindividual variability in the modulation of movement intention by AG stimulation. In line with the secondary role of parietal activity in the 450–950 ms time window after movement, adding this activity to the general linear model yielded no further improvement, which still accounted for 80.9% of behavioral variance.
Peak LRP amplitude correlates with W-time modulation
We did not find significant across-subject correlation between the modulation of W-time by M1 stimulation and its modulation of BP, LRP, or parietal activity after Bonferroni correction (for a potential account, see Computational model). Next we adopted a median-split approach by separating all subjects into an early-W and a late-W group based on their reported W-times in each condition (sham, AG, and M1 stimulation). The reported W-times were significantly different between the two groups in all conditions (all p < 0.005). We then computed the LRP for each group (Fig. 5A,B). The biggest difference between early and late W-time groups resided in the peak amplitude of LRP around button press (assessed as the mean of a 300 ms window centered around movement), which was significantly larger in the early-W than the late-W group under both AG and M1 stimulation conditions (p < 4e-7 and p = 0.014, respectively). In addition, the peak LRP amplitude exhibited a monotonic relationship with the reported W-time across groups and conditions: the larger the LRP amplitude, the earlier the reported W-time (Fig. 5C; Spearman rank correlation, ρ = 1, p < 0.003). Consistently, at the population level, the peak LRP amplitude was significantly larger in both AG and M1 stimulation conditions than sham (p < 0.004 and p < 0.03, respectively), which is in line with W-times being earlier in these two conditions (Fig. 3). Neither BP nor parietal activity showed the same pattern of results.
LRP peak amplitude correlates with W-time. A, Subjects were split into two groups based on their reported W-times in AG anodal and sham conditions, and the mean LRP was calculated for each group. The peak LRP amplitude was assessed in a 300 ms window centered around movement (yellow window). Mean reported W-times for each group are marked as ticks along the x-axis. B, Same as A, for M1 anodal stimulation (stim) compared with sham condition. Sham data are identical to those in A. C, Scatter plot for peak LRP amplitude against W-time for each group of subjects defined by the median splits in A and B. The location of the crosshair and the error bars indicate mean ± SEM across subjects. The W-time is significantly correlated with LRP peak amplitude across subject groups (Spearman rank correlation ρ = 1; p = 0.0028).
Computational model
The above results demonstrate that a temporally extended process lasting seconds around movement onset contributes to the reported time of movement intention. Brain activity as early as >2.5 s before movement and as late as 550 ms after movement contributes independently to the final reported W-time (Fig. 4). A parsimonious account of our results is that the final reported W-time is computed as a weighted average of the visual representation of the revolving clock during this process, and both AG and M1 anodal stimulation shifts this process earlier in time (Fig. 6A). To formalize these intuitions, we constructed a computational model using coupled leaky integrator units to represent activity in M1, AG, and SMA (Fig. 6B,C; for details, see Materials and Methods).
Schematic of the computational model. A, A conceptual account. Our empirical data suggest that W-time is influenced by neural activity occurring ∼2–3 s before, and up to 1 s after, button press. One parsimonious account is that the brain process underlying movement intention lasts for seconds, spanning the period before and after movement, and that tDCS shifts the temporal boundaries and relative intensity of this process across time. Assuming that W-time is generated by a weighted average of the perceived clock time during the intention process, shifting this process earlier in time would lead to an earlier estimate of W-time. The subjective perception of clock time may lag the physical time by ∼80 ms, but this lag is assumed to be constant across tDCS conditions. B, Model structure and all relevant parameters for simulating the ERP results. Three leaky integrator units are simulated, corresponding to SMA/pre-SMA, M1/PMC, and AG, respectively. Their connectivity pattern is shown. W parameters describe the influences between nodes, and their potentiation under tDCS is shown by asterisks (corresponding to the k parameters in the model). Parameters capturing tDCS effects are shown in red. For a detailed model description, see Materials and Methods. C, Schematic time courses of all time-varying parameters in the model.
Simulated M1 and AG activity in the time window [−4 s, 2 s] are shown in Figure 7, B and D, respectively. These activity time courses were subjected to baseline correction (Fig. 7A,C) to compare with empirical LRP and parietal activity, respectively. Our model accounts for the following features of the empirical ERP results: (1) under sham stimulation, the LRP remains close to baseline during [−4 s, −1] and then undergoes a large peak with onsets around t = −1 s and peaks around t = 0 (Fig. 4A, data; Fig. 7A, model); (2) under AG stimulation, the LRP is enhanced in the interval t = [−4 s, −1 s] and the LRP peak around t = 0 is also enhanced (Fig. 4A, data; Fig. 7A, model); (3) under M1 stimulation, the activity of the LRP in the interval t = [−4 s, −1 s] is not enhanced, but the LRP peak around t = 0 is enhanced (Fig. 5B, data; Fig. 7A, model); (4) under sham stimulation, the parietal activity gradually ramps up over the interval t = [−4 s, 0 s], with an inflection point around t = −1 s (Fig. 4C, data; Fig. 7C, model); (5) under AG stimulation, the parietal activity returns to baseline after t = 0 faster than under sham stimulation (Fig. 4C, data; Fig. 7C, model); and (6) under M1 stimulation, the parietal activity does not significantly differ from that under sham stimulation (data not shown; Fig. 7C, model).
Model output for M1 and AG activity. A, B, Simulated xM1(t) time courses before (B) and after (A) baseline correction under different tDCS conditions. Raw simulated xM1(t) activity (shown in B) is used for W-time computation (Fig. 8), while the baseline-corrected activity (shown in A) is used for comparison with the empirical LRP. C, D, Simulated xAG(t) time courses before (D) and after (C) baseline correction. The baseline-corrected activity (C) is used for comparison with the empirical parietal ERF. By supposing that tDCS affects baseline activity and synaptic strengths in the stimulated brain region, the model reproduced the main pattern of ERP results observed in the empirical data. Specifically, the model reproduces the empirical findings that AG stimulation elevates early and peak amplitude of the LRP and causes a faster return to baseline of parietal activity, whereas M1 stimulation elevates only the peak amplitude of the LRP.
An intuitive explanation of how the model captures the effects of tDCS on ERP is as follows. When AG is stimulated, its baseline activity is elevated (Fig. 7D), although the elevated baseline is rendered invisible by baseline correction customarily applied in ERP analysis (Fig. 7C). With an elevated baseline, AG more powerfully activates M1 than it does under the sham condition, which accounts for the elevated LRP in the time window [−4 s, −1 s], as well as the enhanced peak around t = 0 (Fig. 7A,B). Under AG stimulation, AG is also sensitized to the inhibitory input from M1, and this accounts for the faster return to baseline of the parietal activity (Fig. 7C). When M1 is stimulated, its baseline activity is also elevated (Fig. 7B), although the elevated baseline is again rendered invisible by baseline correction (Fig. 7A). Further, M1 stimulation causes M1 to become sensitized to input from SMA, resulting in an elevated peak in the LRP around t = 0, which remains after baseline correction (Fig. 7A).
Using a readout mechanism based on temporal averaging of the visual input, where weights depend on M1 activity (Fig. 8A–C; for details, see Materials and Methods), the model produced W-times of −190 ms in the sham condition, −246 ms in the AG stimulation condition, −263 ms in the M1 stimulation condition, and −187 ms in the SMA stimulation condition, a close match to the empirical data (Fig. 8D). Thus, the model captured the main findings that both AG and M1 stimulation shifted the reported W-time significantly earlier relative to sham, whereas SMA stimulation did not. Intuitively, the success of this model is due to the fact that M1 activity is elevated in the early time period under both AG and M1 stimulation (Fig. 8B), which shifts the temporal weighting function w(t) slightly earlier than in the sham condition (Fig. 8C). This shift lends more weight to early time values in the computation of W-time, resulting in earlier computed W-times.
Simulation of W-time computation and across-subject correlations between tDCS effects on neural activity and W-time. A, Components entering into the W-time computation in the sham condition, normalized to the peak of each time course. xM1(t) is the simulated M1 activity under sham condition. xM1(t) is weighted by a fixed weighting function β(t) to yield the final set of weights, w(t). W-time is computed as a weighted average of perceived clock time in the time interval t = [−2850 ms, 950 ms], using w(t) as the weights (Eq. 8). B, Simulated M1 activity under AG and M1 stimulation and sham conditions, without baseline correction (same as in Fig. 7B). C, w(t) under the same conditions. Due to elevated M1 activity in the early portion of the time window under AG and M1 stimulation, as depicted in B, w(t) distribution is shifted slightly earlier than in the sham condition. D, Simulated W-times from the model closely match the mean reported W-times by the subjects under different tDCS conditions. E, Using a Monte-Carlo simulation of 1000 subjects, the model captures the across-subject correlation between the effects of AG stimulation on early LRP and W-time (r = 0.99), as well as the correlations between the effects of AG stimulation on late parietal activity in both temporal clusters and W-time (r = −0.97). F, The simulated 1000 subjects were split into two groups based on their computed W-times under AG (left) and M1 (right) stimulation; the mean M1 activity was calculated for each group and baseline corrected to compare with the empirical LRP (Fig. 5A,B). Simulated M1 activity under the sham condition is shown for comparison.
We next computed the across-subject correlation between tDCS modulation of W-time and its effects on early M1 activity and late AG activity, using the same method as applied to the empirical data (Fig. 4). The simulation yielded an early LRP–W-time correlation of r = 0.99, and a late parietal–W-time correlation of r = −0.97 for both parietal clusters in the time windows of 0–550 and 450–950 ms (Fig. 8E). Thus, the model captures the correlation patterns observed empirically. The empirical correlations are slightly lower, with an early LRP–W-time correlation of r = 0.71, and a late parietal–W-time correlation of r = −0.78 and −0.75, respectively, for the first and second time windows, presumably due to noise in the recording system and potentially additional neural factors not accounted for in the model.
Last, we investigated whether the model could capture the empirical result of a larger peak LRP amplitude in the early-W group than in late-W group under both AG and M1 stimulation (Fig. 5). The simulated 1000 subjects were split into two groups based on their W-times under AG and M1 stimulation conditions. The simulated M1 activity (after baseline correction, for comparison with the empirical LRP) was averaged within each group and shown in Figure 8F. Since variability was only introduced to tDCS-modulated parameters, sham results were identical across the simulated subjects. These results show that the peak amplitude of simulated LRP is indeed higher in the early-W group than in the late-W group under both AG and M1 stimulation, capturing the patterns seen in the empirical data.
Testing model predictions
Above we have shown that a computational model comprising coupled leaky integrators reproduces the main findings in the empirical data, including tDCS modulation of W-time and ERPs, as well as the correlations between these effects. This model postulates a simplified connectivity pattern among three nodes—SMA/pre-SMA, M1/PMC, and AG—that is inspired by prior experimental findings (see Computational model in Materials and Methods). The key insight from the model is that by incorporating empirically observed neurobiological effects of anodal tDCS, namely, enhanced spontaneous neural excitability and synaptic potentiation (Nitsche et al., 2008), and assuming a simple temporal-averaging computational mechanism of W-time, the model was able to reproduce the full pattern of empirical findings regarding tDCS modulation of W-time and the associated neural changes. While it is impossible to rule out all alternative mechanisms to explain our empirical result showing that early LRP and late parietal activity independently correlate with the reported W-time, an extended temporal-averaging process appears to be the most parsimonious explanation. Last, we wondered whether it is possible to directly observe model predictions about the effects of tDCS in our empirical EEG data.
First, in line with earlier physiological studies (Creutzfeldt et al., 1962; Purpura and McMurtry, 1965), the model assumes that anodal tDCS enhances spontaneous neural excitability, which predicts a shift in the baseline activity of LRP under M1 stimulation and in the baseline activity of parietal ERP under AG stimulation (Fig. 7B,D). This baseline shift is invisible in our empirical ERP results (Figs. 4E, 5B) due to the baseline correction procedure used in EEG signal preprocessing (see Materials and Methods). Baseline correction is typically considered necessary for analyzing slow ERPs due to drift artifacts in DC-EEG recordings (Rockstroh et al., 1989). Nonetheless, could we observe a baseline shift in the EEG signal due to tDCS as predicted by the model? An additional ERP analysis following a procedure similar to that described above, except without baseline correction, suggested that this is indeed the case—there was a baseline shift in the LRP under M1 stimulation and, similarly, in the parietal ERP under AG stimulation, as predicted by the model (Fig. 9B,D). Baseline-corrected ERPs as presented earlier, are reproduced in Figure 9, A and C (with arrows indicating effects that correlate with the change in W-time), to compare them with model output (Fig. 7).
Testing model predictions. A, C, Empirically observed LRP and Parietal ERP as a function of tDCS condition. Arrows indicate the explanatory targets for the model: portions of the LRP and parietal ERP whose modulation by tDCS is correlated with tDCS-induced shifts in W-time (Figs. 3, 4, 5). B, D, Same as A and C, but computed without performing baseline correction. The model predicts that M1 and AG stimulation elevates baseline activity in these regions (Fig. 7B,D); a qualitatively similar pattern is seen in the empirical ERPs without baseline correction. E, Functional connectivity between LRP and parietal activity electrode groups (corresponding to the M1/PMC–AG pathway in the model) under sham, AG anodal, and M1 anodal conditions. Functional connectivity was computed as the Pearson correlation on spontaneous activity extracted from the intertrial interval periods and Fisher z-transformed. F, Functional connectivity between LRP and BP electrode groups (corresponding to the SMA/pre-SMA–M1/PMC pathway in the model) under sham and M1 anodal conditions. In all panels, group-averaged results across 12 subjects are shown.
The second model assumption regarding the effect of tDCS is synaptic potentiation, which has been previously observed in physiological studies (Fritsch et al., 2010). We reasoned that such synaptic potentiation might manifest as changes in functional connectivity between brain regions that could be probed using our empirical EEG data. It is important to stress that the model posits a simplified pattern of connectivity inspired by prior experimental findings, which constituted a minimal set of assumptions required to capture the key empirical findings in this study. Functional connectivity in EEG data, on the other hand, is influenced by many other factors not accounted for in the model, such as input from other brain regions not simulated in the model (for a discussion on model simplifications, see Computational model in Materials and Methods). Nonetheless, we computed functional connectivity between spontaneous EEG signals from different electrode groups using the intertrial interval period (see Testing model predictions in Materials and Methods) and observed trend effects that were consistent with model predictions (Fig. 9E,F). First, the functional connectivity between LRP and parietal electrode groups was reduced under AG stimulation (Fig. 9E), which is consistent with the model postulation of an inhibitory pathway from M1 to AG that would be potentiated under AG stimulation (Fig. 6B). By contrast, LRP–parietal connectivity was enhanced under M1 stimulation (Fig. 9E), which is consistent with the model assumption of an excitatory input from AG to M1 that would be potentiated under M1 stimulation (Fig. 6B). Last, the connectivity between LRP and BP was enhanced under M1 stimulation (Fig. 9F), accordant with the model postulation of an excitatory input from SMA to M1 that is potentiated under M1 stimulation (Fig. 6B). These qualitative observations provide further support for the connectivity pattern and tDCS effects posited in the model, which were inspired and constrained by prior empirical findings.
Discussion
We observed that anodal HD-tDCS over both M1 and AG significantly altered the reported W-time during a self-generated movement task. In both cases, the reported W-time became substantially earlier, moving from −188 ms under sham condition to −249 ms under AG stimulation (p = 0.016), and −258 ms under M1 stimulation (p < 0.009; Fig. 3A). By contrast, SMA stimulation did not significantly alter W-time (p > 0.05), and none of the stimulation conditions significantly altered the reported M-time or the reported S-time.
While we are not aware of any precedent in the literature showing that M1 modulation influenced movement intention, previous studies have shown that patients with lesions in AG report W-times that are much later than those reported by healthy subjects (Sirigu et al., 2004), and direct cortical stimulation of AG can produce conscious intention to move and even illusory movement perception (Desmurget et al., 2009). These results are consistent with our observation that enhancing neural excitability in AG using anodal tDCS moved the reported W-time significantly earlier. Together, these findings support a recent framework differentiating the involvement of AG in the initial formation of movement intention from the role of SMA/pre-SMA in releasing M1 from inhibition more proximal to movement onset (Desmurget and Sirigu, 2009, 2012). Nonetheless, it remains a possibility that our negative finding with SMA/pre-SMA stimulation was due to insufficient current amplitude within HD-tDCS safety parameters given the cortical depth of this region (Fried et al., 1991, 2011). An MRI-guided finite element model (FEM) of the current distribution under HD-tDCS showed that the spatial focality of HD-tDCS is vastly improved over traditional bipolar tDCS and, as expected, current density falls off with increasing cortical depth (Datta et al., 2009, 2012). This model was partially validated using empirical data (Datta et al., 2013; Edwards et al., 2013). In particular, an FEM simulation using a 6 cm ring montage (similar to our 5 cm ring montage) found that with M1 stimulation, current distribution reaches as deep as the ventricles (Datta et al., 2012). Since current conductivity along the medial wall is higher than through brain tissue, due to the higher conductivity of CSF and blood vessels, it is very likely that HD-tDCS current under our stimulation parameters reached SMA/pre-SMA. Nonetheless, future FEM simulation of HD-tDCS over medial brain regions will be helpful for interpreting the current results.
To avoid electrical interference between tDCS and EEG recording, subjects performed the self-generated movement task under EEG recording right after the termination of tDCS current. On average, the experiment was performed within a window of 24–71 min after the termination of tDCS current. A previous study has established that the aftereffect of HD-tDCS delivered at a lower intensity (2 mA, 10 min) than that used in the current study (2 mA, 20 min) lasts for >2 h. Hence, the modulation of W-time by tDCS observed in this study was entirely mediated by its modulation of spontaneous neural activity and synaptic strength. Recent neuroimaging and electrophysiology studies have demonstrated the functional and cognitive relevance of spontaneous brain activity (Raichle, 2006; Boly et al., 2007; Fox et al., 2007; Hesselmann et al., 2008; Wyart and Tallon-Baudry, 2009; Berkes et al., 2011; He, 2013); the present work, using a causal manipulation, adds significantly to this growing literature.
Our EEG recording obtained during task performance revealed the neural mechanisms underlying tDCS modulation of W-time. We observed that AG stimulation enhanced LRP amplitude 1–3 s before movement onset, and caused parietal activity to return to baseline faster after movement (Fig. 4). Strikingly, these neural effects independently contributed to tDCS modulation of W-time, and together accounted for 81% of across-subject variance in the behavioral effect. In addition, both AG and M1 stimulation produced a larger LRP peak amplitude around movement onset, specifically in subjects with the strongest W-time modulation (Fig. 5). Mechanistically, how does tDCS modulation of spontaneous neural activity induce changes in ERPs under task performance, which are time locked to movement onset and baseline corrected? In particular, baseline correction, while a necessary step in EEG signal preprocessing to remove slow artifacts, eliminates any potential difference in baseline neural activity between conditions. This question is addressed by our computational model. Previous studies (Nitsche et al., 2008; Fritsch et al., 2010) have established that anodal tDCS induces depolarization of resting membrane potentials and slightly elevated spontaneous neuronal firing rates, as well as increased synaptic efficacy. We incorporated these mechanisms into the model, which described M1, AG, and SMA as interacting leaky integrator units whose baseline (i.e., spontaneous) activity and pattern of synaptic connectivity are modulated by tDCS (Fig. 6). This model reproduced the empirically observed effects of AG and M1 stimulation on ERPs (Fig. 7). It illustrates how, under recurrent network architecture and with modulation of synaptic strength, tDCS modulation of spontaneous neural activity could translate into changes in ERPs that are time locked to movement onset, after baseline correction. Moreover, the model proposes a potential account for the empirical observation that while M1 stimulation induced few changes in the ERPs, it nevertheless produced a significant earlier shift in the reported W-time. In the model, this is because M1 stimulation induced tonically elevated M1 activity, which influences the computation of W-time (Fig. 8B,C) but is rendered invisible by baseline correction in the ERP analysis (Fig. 7A). Reassuringly, a post hoc ERP analysis without the baseline correction procedure revealed baseline shifts in the ERP due to tDCS stimulation that are consistent with model predictions (Fig. 9B,D).
A simple hypothesis for the genesis of the W-time is that when the underlying neural activity crosses a fixed threshold, movement intention reaches conscious awareness, and that instant constitutes the later reported W-time (Lau et al., 2006; Fried et al., 2011). While the race-to-threshold model is an attractive framework for perceptual decision making (Gold and Shadlen, 2007) and movement execution (Schurger et al., 2012), it cannot explain our (Fig. 4C–F) or previous TMS and psychophysics results (Lau et al., 2007; Banks and Isham, 2009) on conscious movement intention. These results show that brain activity happening hundreds of milliseconds after movement can still influence W-times reported to be much earlier. By contrast, in line with the hypothesis that conscious movement intention is generated by a process that unfolds in parallel with the unconscious initiation of movement (Wegner, 2003; Hallett, 2007; Lau et al., 2007; Desmurget and Sirigu, 2009), our results reveal a temporally extended process contributing to the final reported time of movement intention (illustrated in Fig. 6A). This process starts seconds before movement (as supported by our early LRP result; Fig. 4A,B) and concludes hundreds of milliseconds after movement (as supported by our late parietal result; Fig. 4C,D). The final reported time of movement intention likely depends on ongoing neural computations occurring during this process, perhaps consisting of a weighted average of the visual representation of the revolving clock. In addition to the start and end times of this process, its intensity may also play a role in the final reported W-time (our peak LRP result; Fig. 5), potentially by changing the weight profile of the temporal averaging (Figs. 6A, 8C). It is important to stress that this account does not in any way imply physical backward causality, because although the W-time was reported to be ∼200 ms before movement, the report itself was made much later (2.56–5.12 s after movement in this study).
We captured these qualitative observations in the model, which implements in a quantitative way the above account of how W-time might be computed from a weighted average of the visual input over a time period spanning from ∼2.8 s before to ∼1 s after movement, where the weights depend on M1 activity (assumed to underlie the LRP). With these minimal assumptions, the model reproduced the earlier reported W-times under both AG and M1 stimulation conditions (Fig. 8D). To simulate individual subjects, we introduced variability to the parameters describing the effects of tDCS. The model output captured all main aspects of the empirical data, including across-subject correlations between tDCS modulation of the reported W-time and its effects on the early LRP and late parietal activity (Fig. 8E), as well as higher LRP peak amplitude in the early-W group than the late-W group (Fig. 8F). This model thus illustrates the plausibility of a temporally weighted averaging process contributing to the reported time of movement intention. Mechanistically, such temporal averaging could potentially be performed by trajectory-based processing that is time sensitive (Maass et al., 2002; Rabinovich et al., 2008; He, 2013).
Last, in addition to the parallel, temporally extended processes underlying intention and movement initiation (Fig. 6A), there may be another, third, process that drives both of them. This third process might correspond to activity in high-level prefrontal and posterior cingulate cortices (Soon et al., 2008), or an unconscious “prior intention” process, as previously proposed (Desmurget and Sirigu, 2009). In principle, tDCS in the current study could have differentially modulated the influences of this third process on the intention and movement initiation processes. Although this scenario is unlikely given that we did not directly stimulate these high-level areas and given the relatively crude spatial resolution of tDCS, at present it cannot be entirely ruled out.
In conclusion, we observed that enhancing spontaneous neuronal excitability in both M1 and AG substantially lengthened the interval between the reported time of movement intention and the time of movement execution. Slow brain waves recorded by EEG explained as much as 81% of interindividual variability in this effect, revealing a temporally extended process spanning seconds around movement onset that contributes to the final reported time of movement intention. In addition, a computational model using biologically based assumptions about the effect of tDCS on brain networks accounted for all main aspects of the empirical findings. While the possibility of a temporally extended process underlying movement intention was previously raised (Lau et al., 2007; Matsuhashi and Hallett, 2008), the present data are the first to directly demonstrate such a process. Interestingly, patients with psychogenic tremor report W-times much closer to movement onset than control subjects (Edwards et al., 2011), and have reduced activity in AG (Voon et al., 2010). Hence, our results might suggest potential avenues of treatment for neuropsychiatric illnesses involving a distorted sense of volition. Last, in line with a recent study (Li et al., 2014) showing that slow cortical potentials (SCPs) correlate with conscious visual perception, here we observed that causal manipulations of SCPs influenced the reported time of conscious movement intention, which further extends the link between SCPs and conscious awareness (He and Raichle, 2009).
Footnotes
This research was supported by the Intramural Research Program of the National Institutes of Health/National Institute of Neurological Disorders and Stroke. We thank Irene McClain for medical coverage; Adam Steel for assistance with the Brainsight system; and Qi Li, Alex Baria, and Amy Lin for help with setting up EEG electrodes.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License Creative Commons Attribution 4.0 International, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
- Correspondence should be addressed to Biyu J. He, 10 Center Drive, Building 10, Room B1D728, Bethesda, MD 20892-1065. biyu.jade.he{at}gmail.com
This article is freely available online through the J Neurosci Author Open Choice option.
