Can electric fields explain inter-individual variability in transcranial direct current stimulation of the motor cortex?

The effects of transcranial direct current stimulation (tDCS) on motor cortical excitability are highly variable between individuals. Inter-individual differences in the electric fields generated in the brain by tDCS might play a role in the variability. Here, we explored whether these fields are related to excitability changes following anodal tDCS of the primary motor cortex (M1). Motor evoked potentials (MEPs) were measured in 28 healthy subjects before and after 20 min sham or 1 mA anodal tDCS of right M1 in a double-blind crossover design. The electric fields were individually modelled based on magnetic resonance images. Statistical analysis indicated that the variability in the MEPs could be partly explained by the electric fields, subjects with the weakest and strongest fields tending to produce opposite changes in excitability. To explain the findings, we hypothesized that the likely locus of action was in the hand area of M1, and the effective electric field component was that in the direction normal to the cortical surface. Our results demonstrate that a large part of inter-individual variability in tDCS may be due to differences in the electric fields. If this is the case, electric field dosimetry could be useful for controlling the neuroplastic effects of tDCS.

Scientific RepoRTs | (2019) 9:626 | DOI: 10.1038/s41598-018-37226-x experiments, we applied 1 mA anodal tDCS for 20 min on the right M1. Previously, a similar protocol has been shown to decrease the excitability of the right M1 21 . A similar inhibitory effect of long-duration anodal tDCS has also been shown using 26 min of 1 mA on the left M1 22 . Notably, halving the duration to 13 min enhanced the excitability 22 , consistently with the typical facilitatory effect of 9-13 min long anodal stimulation 23 .
To find which cortical sites are potentially affected by the EF, we decided to use partial least squares (PLS) regression 24,25 , which is an effective method for finding relationships between dependent variables (here: MEP amplitude) and a large number of collinear predictor variables (here: EF in the cortex). Compared to other commonly used approaches for feature extraction from imaging data, such as random field theory 26 , PLS regression was advantageous because we needed not define a region of interest a priori, which would have been arbitrary as we did not know in advance which site in M1 or other regions 27 was affected by tDCS. At the potentially important cortical site, the data were further analysed using linear mixed effects models to investigate the direction and time-dependence of the effects.

Methods
Subjects. Twenty-eight subjects (7 females and 21 males; mean age ± SD = 27 ± 6 years; 26 right and 2 left handed) participated in the experiments. The subjects were the same who participated in our previous study 28 . The subjects were neurologically healthy and had no family history of epilepsy. The handedness was assessed using the Oldfield handedness questionnaire 29 . All subjects gave informed consent before participating in the experiments. The Human Ethics Committee at the National Institute for Physiological Sciences, Okazaki, Japan, approved the experiments. All methods were carried out in accordance with approved institutional guidelines and regulations.
Experimental parameters. The experiment employed a double blind, sham-controlled, crossover design to study the effects of anodal tDCS over the right M1 on the MEPs. The experimental parameters are summarized in Fig. 1.
TDCS (1 mA) was applied using a DC STIMULATOR PLUS (neuroConn GmbH, Ilmenau, Germany) with two conditions. Conditions were separated by a washout period of at least three days. The MRI of each subject were acquired prior to the experiments.
The stimulation (anode) electrode (surface area 5 × 5 cm 2 ) was placed over the hand M1 in the right hemisphere. The cathode (surface area 5 × 5 cm 2 ) was placed over the contralateral orbit. Stimulation lasted for 20 min. In the sham condition, current (1 mA) was applied only for the first 15 s. The fade-in/fade-out time was 10 s in each condition. The subjects were asked to sit on a chair during the experiments, and the experimenter observed that the subject was maintaining rest.
The location of the anode was identified using an individual T1-weighted MRI and a frameless stereotaxic navigation system (Brainsight 2; Rogue Research Inc., Montreal, Canada). The experimenter first identified the "hand knob" structure of the precentral gyrus, the centre of which was projected to the scalp, which is illustrated in Fig. 1B. The centre of the anode was placed at the projected point on the scalp. TMS was also applied at the same point, in the direction perpendicular to the central sulcus, as identified using the navigation system.
As a measure of cortical excitability, MEPs were elicited using a Magstim 200 2 magnetic stimulator (The Magstim Company Ltd, Whitland, UK). At the beginning of each condition, we determined the resting motor threshold (RMT). RMT was defined as the lowest stimulation intensity required for eliciting MEPs of 50 μV peak-to-peak amplitude in five of ten trials in the fully relaxed left abductor pollicis brevis (APB) muscle 30 . The TMS coil was held in place manually, and its position was kept constant using the navigation system.
MEPs of the fully relaxed left APB muscle were recorded before and 0-30 min (with 10 min intervals) after tDCS. During test stimulation, TMS with the intensity of 130% of the RMT was applied 30 times for each time point. The TMS pulse interval was randomly assigned at 6, 7, 8, 9, or 10 s. Recorded electromyograms were amplified (1000x), band-pass filtered (10-1000 Hz), and sampled at 5 kHz. Finally, the mean peak-to-peak MEP amplitude was calculated.
Anatomic models and inter-subject registration. T1-and T2-weighted MRI were segmented into distinct tissue compartments. Brain tissues were segmented using the FreeSurfer image analysis software [31][32][33][34] , and the remaining tissues were segmented using custom methods implemented in MATLAB (The MathWorks Inc., Natick, MA, US). The segmentation process and the tissue conductivities were identical to our previous studies 16,28 . The conductivities were (unit: S/m): grey matter 0.2, white matter 0.14, blood 0.7, compact bone 0.008, spongy bone 0.027, CSF 1.8, dura and muscle 0.16, skin and fat 0.08, and eye 1.5.
FreeSurfer with the default parameter values was used to generate a mapping from the surface of each individual subject's brain to that of the standard brain; details of the procedure have been described earlier 16 . The standard brain was based on Montreal Neurological Institute (MNI) ICBM 2009a nonlinear asymmetric template 35,36 . Electric field modelling. The EFs calculated in each subject were identical to those reported in our previous study 28 . Briefly, the EFs were modelled using the following methods. The electrodes were modelled using a  16 . The electrical sources were a current source (1 mA) and sink (−1 mA) placed inside the rubber pad of the anode and cathode, respectively. The FEM with cubical 0.5 mm × 0.5 mm × 0.5 mm first-order elements was used to determine the electric scalar potential φ from the Laplace-type equation σ φ ∇ ⋅ ∇ = 0. The equation was numerically solved using the geometric multigrid method 37 to the relative residual of 10 −6 .
In each subject, the EF was calculated from φ → = −∇ E at the depth of 1 mm below the grey matter surface. We then calculated the absolute value (E abs ) and the normal component of the EF where → n is the inner normal vector of the cortical surface).
In order to compare the EFs from different subjects, E abs and E n were mapped to the MNI brain 16 . The final EFs were represented on a triangular surface mesh of the MNI brain, which consisted of 149319 vertices in the right hemisphere and 148076 vertices in the left hemisphere. Data analysis. MATLAB (version 2017b, The MathWorks, Inc.) was used for all statistical tests. The significance level was P < 0.05. Outliers were detected using Grubbs' test.
Effects of Time, Session and EF on the MEP. We first analysed the experimental results without considering the EF, as would conventionally be done in tDCS studies. A linear mixed effects model was used to study the effects on the MEP normalized to the baseline. As fixed effects, we entered the effects of Time (t = 0, 10, 20, and 30 min after stimulation, denoted t0-t30), Session (real tDCS and sham) and their interaction effect. To account for the effects of the baseline MEP (MEP base ) on the normalized MEPs, the fixed effects of MEP base , MEP base × Time, MEP base × Session, and MEP base × Time × Session were also included in the model. By-subject intercepts were treated as random effects. Maximum likelihood estimation was used for determining the linear mixed effects model parameters.
To study the effect of the EF on the normalized MEPs, the absolute value and normal component of the EF were calculated at an observation point, → r 0 (to be defined later). We then added additional fixed effect terms for EF, Time × EF, Session × EF, and Time × Session × EF to the linear mixed effects model. The likelihood ratio test was used to test whether adding the EF terms improved the model significantly.
Using the linear mixed effects models, we studied the following questions: (1) was the mean value of the normalized MEP different from the baseline at any time point, (2) did the mean value depend on Session, (3) was the slope for MEP base nonzero at any time point and (4) did it depend on Session, (5) was the slope for EF nonzero at any time point and (6) did it depend on Session, and (7) did any of the mean values or slopes depend on the time point? The coefTest function of MATLAB, which uses F-tests, was used for calculating the P-values. The Satterthwaite approximation was used for estimating the degrees of freedom. For visualization of the linear mixed effect models, we used a similar approach to calculate the P-values for each mean value, slope, and the difference in slopes between sham and real tDCS.
To test the effects of Session and EF on MEP base , we used paired two-tailed t-tests and/or Pearson correlation coefficients.
As a measure of the overall change in the cortical excitability, we calculated the mean MEP amplitude normalized to the baseline over post-stimulation time points (t0-t30) Before this calculation, we verified that there were no significant differences between the post-stimulation time points.
We have previously found in the same subjects that the RMT correlated with the EF strengths in hand M1, and therefore had potential as a simple estimate of the EF strength 28 . To investigate whether the RMT was useful for estimating the normalized MEPs, we repeated the analysis replacing the EF with the RMT.
Estimation of important brain regions using PLS regression. We used PLS regression in MATLAB to study whether the measured MEPs could be explained using the calculated EFs, and, if yes, which brain regions were important for the prediction.
The input data to the PLS regression model were the following. The predictor variables, matrix X (28 × n), were the highest EF values (E abs or E n ) on the right hemisphere. The n vertices with the highest EF values were selected by calculating the (100 − r E )th percentile of the mean value of E abs or |E n | over all 149319 vertices. To test the robustness of the approach, r E was varied from 1 to 10%. The dependent variable, vector Y (28 × 1), was the mean normalized MEP. The columns of X were scaled by dividing them by their sample standard deviations and centred by subtracting their sample mean.
In the initial analysis, the number of PLS components (not to be confused with EF components) was varied from one to five, and the goodness of fit was measured in terms of R 2 (multiple correlation coefficient) and Q 2 (cross-validated R 2 ). R 2 and Q 2 are the upper and lower bounds, respectively, of how well the model explains the data and predicts new observations 25 . To calculate Q 2 , we used 10-fold cross validation with 1000 Monte-Carlo repetitions. PLS component i was defined to be predictively significant if . The analysis was also repeated for the sham MEP data. In this case, the EF should be unrelated to the MEP, and thus, no PLS components should have predictive significance.
After finding the number of significant components, we calculated the variable importance for the projection (VIP) to identify which brain regions were important for predicting MEPs from the EFs. Based on the important variables of PLS regression, we selected a single observation point, → r 0 , in an anatomically relevant location to interpret the effect of the EF on the MEP amplitude using linear mixed effects models.

Results
None of the participants reported side effects. Despite no significant differences between sham and real tDCS at the group level, these results indicated that stimulation had some effect, as the normalized MEPs of sham and real tDCS did not correlate within individuals, and the difference between the MEPs of sham and real tDCS after stimulation was larger than that expected from the baseline measurements. This could mean that the size and/or direction of the response differed between individuals. Were these differences due to chance or due to some systematic factor, such as the EF?
Calculated EFs and PLS regression. The absolute value and normal component of the EF were calculated individually in 27 subjects and registered to the standard brain. Figure 3 shows the average EFs in the right hemisphere.
PLS regression analysis with either E abs or E n as the predictor gave up to one predictively significant PLS component, depending on the percentage of the highest EF values used for the analysis (  Fig. 4A. PLS score plots indicated no violations of homogeneity or curvature of the data. Normal probability plots were used to verify the normality of residuals, and no outliers were detected in residual plots either by visual inspection or by Grubbs' test. Next, we investigated which brain regions were important for predicting the mean normalized MEPs using the VIP as the measure. The regions with high VIP in Fig. 4B are candidates for the site of action where the EF has an effect on the MEPs. For the subsequent analysis, we selected the point with the maximum VIP as an observation point, → r 0 , as shown in Fig. 4B. The same point was the global maximum for r E ≤ 5% and a local maximum for larger r E . In MNI coordinates, → = − r (42,13,66) 0 . Based on the probabilistic cytoarchitectonic atlas of FreeSurfer, → r 0 is located at the border between Brodmann areas 6 and 4a. The choice of → r 0 from among several candidates is motivated by the fact that → r 0 is close to the TMS hotspot of the ABP muscle projected to the cortex, (−41 ± 4, −16 ± 4, 60 ± 4), measured in the left hemisphere 39 .
Effects of the EF. Next, we investigated how the EFs at → r 0 affect the normalized MEPs and their time course.
r 0 89), indicating that the EF was approximately normal to the cortical surface at → r 0 . For this reason, we chose to focus on → E r ( ) n 0 in the following analysis. First, we tested the baseline effects of EF. The Pearson correlation coefficients between → E r ( ) n 0 and MEP base were = . r 0 15 ( = . P 0 4) and = . r 0 26 ( = .  Figure 5A visualizes the fitted linear mixed effects model. Grubb's test for the model residuals detected one outlier data point (exceptionally low normalized MEP at t0 in the subject with the lowest EF). Inclusion or exclusion of the outlier did not change the conclusions, and therefore we have included it in the analysis.
F-tests showed that, similarly to the model without the EF, the mean values of normalized MEPs changed significantly from the baseline [ . = . P 0 007]. Investigation of the slopes for → E r ( ) n 0 showed that real tDCS tended to have more negative slopes than that at baseline or those of sham (Fig. 5A). The slopes for sham were not significantly different from the baseline at any time point.
Time did not have significant effects on any of the intercepts or slopes [ . = .   Figure 5B visualizes the linear regression model as well as the individual data. It can be seen that the individuals with the largest EFs showed decreased MEPs compared to sham, whereas the subjects with the lowest EFs showed either absence of effect or increased MEP compared to sham.
Effects of the RMT. Similarly to our previous study in the same subjects 28 , both →   r 0 04, = . P 0 8), indicating that their correlation was due to the EF. For sham, the correlation coefficient between the RMT and the mean normalized MEP was = .
r 0 25 ( = . P 0 2), and the partial correlations without MEP base and the EF were = .

Discussion
We studied the effect of 20 min 1 mA anodal tDCS on the excitability of the right M1, and modelled the EF in each individual subject. The initial findings suggested no significant differences between sham and real tDCS in the group-level MEPs. However, the responses to sham and real tDCS differed at the individual level. Using regression analysis, we showed that individual differences in the EF could partly explain the variability in responses.  If the EF has an effect on the individual MEPs, from which brain regions does the effect originate from? We explored this question using PLS regression to find relationships between the absolute value or normal component of the EF and the MEPs. The results showed that it was possible to predict a part of the effect using linear regression involving the EF normal component. Furthermore, the EF values in the hand area were important for the prediction. While these results are not a direct proof of the exact origin of the effects, the hypothesis that the possible effect originated from the hand M1 is attractive due to several reasons. Firstly, the sites with high importance were located in the lateral part of the hand knob at the anterior bank of the central sulcus [MNI coordinate (±42, −13,66)], which is at the border between the primary motor cortex (BA4a) and premotor areas (BA6). This is in the immediate vicinity of the TMS hotspot of the studied APB muscle 39 , which could explain why the EFs at this site are important for predicting the changes in MEPs. Furthermore, the direction of the EF was approximately normal to the cortical surface. In previous in vitro studies, EFs applied in the normal direction have been shown to produce long lasting after-effects 17,18 . Therefore, a plausible hypothesis for explaining our findings is that the after-effects of tDCS are mediated by the normal component of the EF at or near the TMS hotspot. However, this needs to be confirmed in additional studies.
Linear mixed effects models were used to analyse how E n in hand M1 affected the MEPs. The results showed that tDCS, but not sham, changed the slope between the EF and MEP amplitude, indicating that subjects with low and high EFs responded differently to stimulation. TDCS changed the slope to a more negative direction, i.e., subjects with a stronger E n exhibited a larger decrease (or smaller increase) in the MEP compared to sham or baseline than subjects with a weaker E n . A possible explanation for the negative effect of the EF is that, as shown in previous studies, long duration (≥20 min) 1 mA anodal tDCS may decrease the motor cortical excitability 21,22 . The negative slope is consistent with these findings; the inhibitory effect is stronger in individuals with a stronger EF and absent in subjects experiencing weak fields in M1. However, the effects of tDCS are known to be non-linearly dependent on the current intensity [40][41][42] and stimulation duration 22 . Therefore, our findings should not be extrapolated to other experimental conditions.
As far as we know, our study is the first to combine individual EF modelling with tDCS experiments. However, a few recent studies have used EF models to design and analyse the results of tDCS experiments. The effects of different EF components were studied previously by Rawji et al. 43 , who investigated the excitability changes using two electrode montages that produced EFs dominantly in the posterior-anterior (PA) or medial-lateral (ML) directions (left M1, 1 mA, 10 min, FDI muscle, N = 15). They found that a larger EF in the PA direction decreased the excitability, whereas the EF in the ML direction did not 43 . Our experiments are not directly comparable due to different target muscle, duration, and electrode montage. However, common to both studies was that the effective EF component might have been the normal component, E n , which depends on the representation of the target muscle in the curved surface of the hand area. In another study, Fischer et al. 27 showed that a multifocal tDCS montage that produced a weaker EF in the left M1 resulted in a greater increase in excitability than conventional tDCS (2 mA, 10 min, FDI muscle, N = 15). Consequently, they argued that the effects of tDCS on the motor cortical excitability likely originated from regions outside M1 27 . Although obtained using different stimulation parameters, our results and those of Rawji et al. 43 show that the EF strength-response relationship can also be negative, and thus, the findings of Fischer et al. 27 could also be explained by a local effect of EF in M1.
Existing tDCS protocols that have been used since the early 2000s typically apply the same input current to all subjects 6,23,44 . This approach may be problematic, as our results suggested that the subjects with the lowest and highest EF strengths may respond oppositely to the same input current. Therefore, the findings at the group level may become weak or not significant. Indeed, in our experiment, no significant group-level differences compared to sham or baseline would have been found without considering the inter-individual differences in the EF. For comparison, several previous studies have also reported small group-level responses but high inter-individual variability [7][8][9] . If our results generalize to other stimulation parameters, EF models could be used to select the stimulation current individually, which could reduce variability.
Separately from the effect of the EF, we found that MEP base had a significant effect on the normalized MEPs. The effect was similar for both sham and real tDCS: subjects with a higher MEP base tended to decrease the MEP, and subjects with a lower baseline tended to increase the MEP. Wiethoff et al. 7 have also reported a similar negative effect. We believe that the effect directly follows from normalization and is unrelated to tDCS: If the MEP were a random variable, the conditional expectation of the normalized MEP would be a decreasing function of MEP base .
Despite the correlation between the RMT and the EF 28 , we found that the RMT had relatively weak effects on the after-effects of tDCS. The correlation coefficient between the RMT and mean normalized MEP was = .
r 0 43, which was only marginally better than that obtained for sham. The correlation disappeared when the confounding effect of the EF was excluded, indicating that the RMT and and the post-stimulation MEPs were linked solely through the EF. Interestingly, a correlation of similar strength but with the opposite sign ( = − . r 0 47) has been reported previously for 15 min 1 mA anodal tDCS 42 .
Due to the exploratory nature, our study has several limitations that should be considered in future studies. Firstly, we only considered a single stimulation current (directly proportional to the EF) in each subject. Using additional EFs in each subject could have been used to confirm the regression slopes and take into account inter-individual differences in the sensitivity to EF in linear mixed effects models. If multiple currents are used in future studies, non-linearities [40][41][42] should also be considered when selecting the current values. Secondly, repeating the measurements for each condition would have allowed incorporating intra-individual variability in the statistical models. Notably, previous studies have shown relatively low intra-individual variability for anodal TDCS 11,45 . Thirdly, only one electrode configuration was used in each subject. More than one electrode configuration would have increased the amount of input data to the PLS regression model and would have possibly improved the prediction of the affected regions. Fourthly, we unexpectedly found a significant facilitatory effect Scientific RepoRTs | (2019) 9:626 | DOI:10.1038/s41598-018-37226-x of sham stimulation. We are unsure of the reasons, as other recent sham-controlled studies have not shown any significant effects of sham 11,42,43 . The effect of sham persisted at least 30 min after stimulation, which highlights the importance of sham control in future studies. Future studies can also select the observation point or region of interest a priori. The observation point found in this study is useful for the APB muscle of the left hand. For other targets, the TMS hotspot projected on the cortex may be a reasonable choice. Our results can also be used to estimate the required sample size for future studies. Based on simple linear regression between the EF and normalized MEP, a reasonable estimate for the coefficient of determination is approximately 0.35, which would require at least 20 data points (subjects) to ensure a statistical power of 80%.
In conclusion, this exploratory study showed that individually calculated EFs were related to inter-individual differences in the responses to tDCS. A potential hypothesis for explaining our findings is that the individual effects of tDCS are mediated by the normal component of the EF in the hand area of M1, at or close to the TMS hotspot. If the effect is confirmed, EF modelling could be the key for reducing inter-individual variability in tDCS.

Data Availability
The datasets generated during the current study are available from the corresponding author on reasonable request.