Machine learning and individual variability in electric field characteristics predict tDCS treatment response.

BACKGROUND: Transcranial direct current stimulation (tDCS) is widely investigated as a therapeutic tool to enhance cognitive function in older adults with and without neurodegenerative disease. Prior research demonstrates that electric current delivery to the brain can vary significantly across individuals. Quantification of this variability could enable person-specific optimization of tDCS outcomes. This pilot study used machine learning and MRI-derived electric field models to predict working memory improvements as a proof of concept for precision cognitive intervention.
METHODS: Fourteen healthy older adults received 20 minutes of 2 mA tDCS stimulation (F3/F4) during a two-week cognitive training intervention. Participants performed an N-back working memory task pre-/post-intervention. MRI-derived current models were passed through a linear Support Vector Machine (SVM) learning algorithm to characterize crucial tDCS current components (intensity and direction) that induced working memory improvements in tDCS responders versus non-responders. MAIN
RESULTS: SVM models of tDCS current components had 86% overall accuracy in classifying treatment responders vs. non-responders, with current intensity producing the best overall model differentiating changes in working memory performance. Median current intensity and direction in brain regions near the electrodes were positively related to intervention responses (r=0.811,p<0.001 and r=0.774,p=0.001).
CONCLUSIONS: This study provides the first evidence that pattern recognition analyses of MRI-derived tDCS current models can provide individual prognostic classification of tDCS treatment response with 86% accuracy. Individual differences in current intensity and direction play important roles in determining treatment response to tDCS. These findings provide important insights into mechanisms of tDCS response as well as proof of concept for future precision dosing models of tDCS intervention.

RCT Entities:   Population Interventions Outcomes

Introduction

For almost two decades, transcranial direct current stimulation (tDCS) has been instrumental in advancing the knowledge of human brain function by altering neural activity in the brain [1-3]. tDCS has been suggested to alter resting membrane potentials (i.e., modulate or attenuate), influencing neuronal excitability across the local field of stimulation [1-5]. With this, tDCS has shown great promise as a therapeutic intervention in various neurological and psychiatric disorders [6-9]. While the number of tDCS applications within the literature has grown exponentially [1-5], the optimal dosing parameters (e.g., applied current, electrode placement, etc.) that underlie the positive effect of tDCS remain unclear.

Prior research shows that increasing or decreasing the intensity of applied stimulation to the scalp results in corresponding changes in the electric field within the brain [10]. In vitro studies have shown that the intensity component of tDCS current can modulate cortical excitability [11,12]. Experimental and theoretical studies have indicated that tDCS-related electric field intensity is essential for altering neuron resting membrane potentials and modification of synaptic strength (i.e., LTP/LTD) [13-17]. Studies show that increased applied current intensity is associated with increased amplitude of motor evoked potentials (MEPs), suggesting neuronal sensitization of the motor cortex [8,18-20]. Therefore, varying levels of applied current intensity in tDCS may lead to changes in behavioral outcomes.

Furthermore, electrode placements during tDCS can greatly affect the distribution and direction of electric current throughout the brain [21,22]. Whole-cell recordings have demonstrated that electric field orientation is essential to the likelihood of neuronal firing [23]. Neuronal bodies in parallel with the direction of applied electric fields are more susceptible to stimulation responses. Human studies have also highlighted the importance of tDCS current direction for modulating cortical excitability [24,25]. Rawji et al. (2018) evaluated individual effects of tDCS montages that produced electric fields oriented orthogonal or parallel to the motor cortex on the modulation of MEPs. The orthogonal montage was observed to have greater current flow normal to the cortical surface (i.e., current flow in parallel with the dendritic axis of cortical neurons). These authors reported significant alterations in MEPs with this montage compared to sham [25]. These data suggest that the direction of current flow in the tDCS electric field may be strongly correlated with behavioral outcomes of tDCS.

Conventional tDCS typically employs a fixed applied current (e.g., 2 mA) and electrode placement (e.g., F3/F4) across participants [26]. However, the orientation and intensity of the generated electric field within cortical tissue can be dramatically altered by inter-individual anatomical differences. For instance, brain atrophy can reduce the level of current reaching the brain due to an increase in current shunting within cerebrospinal fluid (CSF) [27,28]. Individual skull thickness and subcutaneous adipose tissue volume can also alter voltage delivery due to differences in tissue conductivity [29-32].

MRI-derived finite element models (FEM) can be used to estimate individualized electric field induced by tDCS. Advancements in tissue segmentation tools and automated modeling pipelines [33,34] have enabled more efficient generation of large and complex FEM that would normally require extensive computing power and time. These models have recently been compared to experimental results obtained via a novel in-vivo magnetic resonance electrical impedance tomography technique [35,36] and intracranial recordings [22,37]. While the experimental results showed a strong correlation with computational model outcomes, a large variation across individuals was observed. Therefore, investigating the nuance of electrical distribution in individualized models may provide more insight into inter-individual variability seen with tDCS. However, the size (i.e., millions of voxels across multiple dimensions) and complexity of generated electric fields has made interpretation of the essential current characteristics (i.e., current direction, intensity, etc.) challenging.

Few studies have attempted to systematically investigate these estimates of electric field distributions as a predictor of responses to tDCS [38-40]. All three prior studies reported increases in current intensity associated with increases in the target behavioral response (e.g., self-reported and physiological measures) [38-40]. Antonenko et al. also reported a positive relationship between the current direction normal to the cortical surface and sensorimotor network strength [40]. These studies employed univariate approaches that treat each voxel or region of interest within the brain as an independent predictor of treatment response. At present, no studies have employed multivariate approaches to investigate patterns within the current distribution as a predictor of treatment response.

Supervised machine learning methods (e.g., support vector machines; SVM) constitute a novel approach in neuroimaging to investigate large and complex datasets [41-44]. SVM uses Mercer’s Theorem [45], which allows the representation of high dimensional feature space in a low-dimensional Gram matrix (Equation (4)) – also known as the “kernel trick” [46]. SVM performs multivariate analyses across many voxels to classify patterns of information [41] that can be used to identify individual contributions of current intensity and direction towards behavioral responses. Multivariate classifiers are iteratively trained to search for patterns within the data that best predict a specified prognostic label, such as behavioral response. This is usually achieved within a cross-validation procedure, withholding a different partition of data for each iteration of training. This is a standard approach within the statistical literature and a widely used technique to provide unbiased generalizability to new data samples [42-44,47]. One study used machine-learning of clinical data to predict treatment responses [48]; however, there are currently no studies that utilize machine learning on FEM to investigate the critical components of dosing parameters in tDCS.

Previous work from our research group demonstrated working memory improvements in older adults following tDCS paired with cognitive training [49]. The objective of this study was to apply machine learning and FEM in the same dataset and identify the central determinants of treatment response. Specifically, the present study applied an SVM machine-learning algorithm to investigate the contributions of current intensity and direction, as well as their interaction, for predicting working memory improvements in older adults. Our primary hypothesis is that SVM applied to individualized tDCS current models are capable of classifying tDCS responders and non-responders above chance (i.e., area under the curve > 0.5). In addition, we hypothesized that the interaction of direction and intensity is the most essential dosing feature for predicting behavioral response. These data will provide critical insight to inform tDCS mechanism theory and provide a potential foundation for methods to increase the effectiveness of tDCS applications.

Methods

Structural imaging and behavioral data were sourced from a phase-II pilot clinical trial that employed a randomized, triple-blinded (assessor, interventionist, participant) design (NCT02137122). This approach enabled examination of the combined effects of tDCS with cognitive training on working memory function in healthy older adults [49].

Participants

Fourteen healthy older adults receiving active-tDCS stimulation were selected for further analysis by the current study [mean (sd) age = 73.57 (7.84), mean MoCA = 27.85 (1.79), 7F:7M]. All participants were screened for eligibility based on study inclusion criteria detailed in the prior manuscript [47]. The study protocol was in accordance with the Declaration of Helsinki and approved by the University of Florida’s Institutional Review Board. Informed written consent was obtained from participants prior to study procedures.

tDCS protocol and application

Conventional 1 × 1 tDCS (Soterix Medical, tDCS-CT for clinical trials) was applied using two 5 × 7cm2 pad electrodes presoaked with 2 mL of 0.9% NaCl and 4 mL added per side (10 mL total per sponge) at F3 (cathode) and F4 (anode) location. Participants underwent head measurements using the International 10–20 system to locate F3–F4 locations at each session. Participants were stimulated at 2 mA intensity for 20 minutes with a 30-second current ramp up and down, with a total of 10 stimulation sessions over 14 days. Each session included 40-minute computerized cognitive training for working memory with stimulation delivered during the first 20 minutes. Details of the computerized training tasks are described in the previous publication [49]. No significant effects of unblinding or differences in sensation were found for active vs. sham participants in the parent study [49].

Behavioral tasks

Participants’ working memory was assessed with an in-scanner N-back task only given at baseline and post-intervention. The task paradigm for each run consisted of four blocks of two-back and four blocks of zero-back presented in a randomized order with 20 seconds of rest between blocks. During the two-back, participants viewed uppercase letters, one at a time. A screen with a central crosshair (+) was presented during the inter-trial interval (Fig. 1). The stimuli appeared for 1 second, followed by a crosshair for 3 seconds, providing a 4-second window to make a response. Details of the N-back task procedure are outlined in the prior paper [49]. Participants performed practice on the N-back task (two- and zero-back) outside of the scanner to ensure understanding of the task at both baseline and post-intervention visits. Two-back performance change (i.e., pre-/post-intervention) was analyzed as a composite percent improvement score for accuracy and reaction time .

Fig. 1.

An example of the stimulus presentation in the two-back variant of the N-back task.

Imaging sequences and parameters

Structural T1-weighted MRI scans were obtained using a 32-channel, receive-only head coil from a 3-T Siemens MAGNETOM Prisma MRI scanner. MPRAGE sequence parameters included: repetition time (TR) = 1800 ms; echo time (TE) = 2.26 ms; flip angle = 8°; field of view (FOV) = 256 × 256 × 176 mm; voxel size = 1 mm3.

Computational model construction

Individual T1-weighted images were converted from DICOM to NIfTI using dcm2niix [50] and resampled with the FreeSurfer v6.0.0 image analysis suite (http://surfer.nmr.mgh.harvard.edu/) into a 256 mm3 field of view (RAS orientation), 1 mm3 voxel size. The computational models of current density were computed using the Realistic vOlumetric-Approach to Simulate Transcranial Electric Stimulation (ROAST; https://www.parralab.org/roast/) toolbox [33] with parallel processing on a high performance cluster with 50 CPU cores and 175 GB of RAM provided by the Research Computing at the University of Florida (HiPerGator). The resampled T1 images (256 × 256 × 256, 1 mm3) were individually processed in parallel using ROAST. The segmentation process was carried out in FreeSurfer to classify tissue types into gray and white matter. FreeSurfer segmentations were visually inspected and manually corrected for errors before reprocessing through FreeSurfer – a procedure that has been validated against manual segmentation [51] and histological measures [52]. Segmentations from FreeSurfer were then combined with segmented CSF, bone, skin, and air from ROAST (See Supplemental Fig. 1). Combined segmented tissues were visually inspected for unassigned voxels to ensure every voxel within the head volume was assigned to one of the six tissue types by overlaying individual segmented volumes with their respective T1 images. Default isotropic conductivity values (gray matter: 0.276 S/m; white matter: 0.126 S/m; CSF: 1.65 S/m; bone: 0.01 S/m; and skin: 0.465 S/m) were assigned to each tissue in ROAST [33]. Volumetric meshes of each tissue type were generated using iso2mesh [53]. Boundary conditions were assigned within ROAST and a finite element solver, getDP [54], was used to compute voltage solutions to the Laplace equation, in the meshed models, where V is the electric potential within the volume and σ is the tissue conductivity. Additional MATLAB routines to compute current density (J) from electric field (E) and tissue conductivity (σ) were added in accordance with Ohm’s law: J = σE. Current density (J) is a useful metric for determining the dosage of current (A/m2) induced in the brain from electrical stimulation. Current density also represents a unit of current directly adjustable through alteration of applied stimulation intensity (e.g., 2 mA vs. 2.3 mA) for the purposes of individualized dosing calculations.

The intensity of current at each voxel was calculated with the function:

Direction of current density in each coordinate plane was separated from intensity by deriving the zenith angle, θ, between and the z-axis, and the azimuthal angle, φ, between the projection of onto the xy-plane, and the x-axis:

Electrode placements

Conventional electrodes (5 × 7 cm2) were constructed from 3D captured electrode models as reported in our previous publication [26] and added to the segmented models in ROAST. Individual electrode placement variability per session for each participant can be found in Supplemental Fig. 3. Default conductivity values of 5.9 × 107 S/m and 0.3 S/m were assigned to the electrode and gel layers, respectively [55]. A +2 mA voltage boundary condition was assigned to the anode electrode at the F4 location, while a −2 mA voltage boundary condition was assigned to the cathode electrode at the F3 location (Fig. 2). Precise electrode locations were recorded for each session with a 3D scanner mounted to an iPad. 3D images were converted to a surface mesh via TechMed3D’s MSoft software [TechMed 3D, Quebec City]. More details of this procedure can be found in a previous publication [26]. Electrode information from each session was added to individual segmented tissue volumes prior to creating volume meshes and solving finite element models, totaling 10 unique head models per individual. Individual, as well as grouped, average and variance maps of current density per electrode placement can be found in Supplemental Figs. 5–8. Electrode displacement was computed as the 3D Euclidean distance between modeled electrodes and the ideal location of F3–F4 derived from ROAST [55].

Fig. 2.

A) The F3-F4 electrode montage with the F3 electrode as the cathode (blue) and the F4 electrode as the anode (red). B) A representative image of actual anode (red) and cathode (blue) placements for a single subject across ten sessions. The mean displacement for this participant of the anode and cathode was 1.9 cm (−0.2 St. Dev.) and 1.7 cm (−0.45 St. Dev.), respectively.

Supervised machine-learning

Current density values were computed in native space then masked using individual participants’ white and gray matter voxels to restrict current values within brain region only and reduce the number of features submitted into the classifier. Masked values were then transformed with the University of Florida Aging Brain-587 tissue probability map [28] into MNI space using SPM’s normalise function [56]. Median values of before and after the transformation were computed as a quality check for transformation errors (r = 0:998). Participant classes were determined by separating participants into binary groups based on pre-/post-intervention performance changes on the two-back working memory task above or below the median. The tDCS responder class (n = 7) had an average performance increase of 22%, whereas, the tDCS non-responder class (n = 7) was found to have an average 9% increase in two-back performance (F(1) = 21:02, p < 0.001). No significant differences in age (p = 0:89), sex (p = 0:66), education (p = 0:53), cortical atrophy (p = 0:22), nor any normalized tissue volume (See Supplemental Figs. 1–2) between these groups were observed.

Due to the high dimensionality of MRI data, feature selection was performed on the training data to further reduce the number of trained features. One popular method of feature selection is to filter the features via voxel wise t-tests between classes to select current elements with a significant group-level difference (p < 0:01) as features for the subsequent prediction step [42,57,58]. To classify responders from non-responders, we used Support Vector Machine (SVM), a machine learning algorithm to search for the optimal hyperplane that separates two classes with maximal margin under the assumption of independently and identically distributed (iid) data [59], which is satisfied in this study. Specifically, LIBSVM [41] was used to optimize the objective function:

where C ≥ 0 is a penalty parameter on the training error. A linear kernel K was generated with the function:

Model performance was evaluated across ten permutations of two-level nested stratified cross-validation [60-62]. To elaborate, we began by splitting the data into K-folds (in this case K = 7) and performed an outer cross-validation loop consisting of K iterations. In each iteration, K-1 folds were used as training data with a single fold left out as test data. A second, inner cross-validation loop was then performed on the training data, providing us with optimal hyper-parameter C. Following training, predictions of held out test data, x, were performed with the decision function:

After all K iterations in the outer cross-validation loop were performed, predicted labels of all subjects were compared against ground truth labels to calculate performance metrics. A receiver operating characteristic (ROC) curve of true positive rate against false positive rate was plotted to demonstrate the separability of classes within each model by calculating the area under the ROC curve (AUC). To assess the stability of feature weights, the standard deviation of feature weights across folds were plotted (see Supplemental Fig. 9). For feature weight generation and deployment, a final model was trained on all fourteen current density maps to derive overall classification weights, w (i.e., the model parameters learned by the optimization function during the training phase, cf. Equation (3)). These weights can be applied to independent data from a new subject to predict their tDCS response classification associated with specific observed J-map features in test data. The feature weights at each voxel, representing the relative contribution of each voxel to the classification, were separated by positive and negative weights that predict responders and non-responders, respectively [63]. Positive and negative weights were divided by their respective sum of weights to compute the percent contribution of each voxel toward either positive or negative predictions. To demonstrate specific brain regions that predict working memory improvements, regions of interest (ROIs) were defined using the Harvard-Oxford atlas [64] and ranked based on their average voxel percent contribution.

Statistical analyses

SPSS Statistics 25 (https://www.ibm.com/analytics/spss-statistics-software) and the Statistics and Machine Learning toolbox in MATLAB 2019b [65] were used to carry out statistical analyses. One-way ANOVA was used to assess mean difference in model performance between the three data types (direction, intensity, direction x intensity). A secondary set of analyses was aimed at determining the characteristics of current within the voxels determined to be essential for treatment response. Within these regions, Pearson’s correlation analyses were used to assess the relationship of behavioral response with field characteristics. Linear regression analyses were used to fit lines of least square residuals. Hedge’s g was computed to define effect sizes of mean differences, corrected for small sample bias [66]. Since all fourteen participants in our study were individuals with no familial relationship and each participant’s data were collected under the same condition, these data points met the statistical assumptions of independently and identically distributed (iid) data. To determine the normality of each variable, we tested the null hypothesis that each data vector comes from a standard normal distribution, against the alternative that it does not come from such a distribution, using the one-sample Kolmogorov-Smirnov test [67,68]. The null hypothesis of a normal distribution was not rejected by the Kolmogorov-Smirnov test for all variables analyzed.

To quantify the required input current for converting non-responders into the responder range, the formula: where x represents a non-responder mean current value, represents the target current values (i.e., average values computed in the responder group), I represents the original injected current (e.g., 2 mA), and represents the percent difference. The optimal electrode displacement was defined as the mean responder displacement minus the non-responder displacement.

Results

Machine learning predictions of tDCS intervention efficacy

Computational models of current intensity, current direction, and their interaction, all demonstrated above chance performance for classifying treatment responders (i.e., AUC>50%). The AUC revealed that the probability of current intensity model to rank a randomly chosen responder higher than a randomly chosen non-responder was 80.6% (Fig. 3A) [69]. Classification of the direction alone and combined models of direction and intensity produced AUCs of 77.6% and 74.9%, respectively. A summary of model performances can be found in Table 1. Computational models of current intensity alone marginally outperformed current direction and the combination of current direction with intensity in the classification problem (Fig. 3B). However, a one-way ANOVA between the AUC of each current characteristic (F(2) = 1:31, p = 0:288) revealed a non-significant difference between the three variables. The support vector machine classification of all three models correctly differentiated tDCS responders from non-responders with averaged overall accuracy of 86.43%. The 95% confidence interval of classification accuracy for these models was between [CI: 85.03–87.83%].

Fig. 3.

A) Areas under the ROC curve across ten iterations of three data types: direction, intensity, and combined. Mean AUC plotted as bars. No significant difference was observed. B) The mean ROC curve of each model across ten iterations with shaded area conveying 95% confidence interval.

Table 1.

A summary table of SVM model performance metrics using ROAST tissue conductivity values [33]: Average Test Accuracy (ACC) with [95% Confidence Interval], Area under the ROC curve (AUC), F1 Score (F1), and Matthews Correlation Coefficient (MCC) for intensity, direction, and direction x intensity of current density.

Intensity				Direction				Direction x Intensity
ACC	AUC	F1	MCC	ACC	AUC	F1	MCC	ACC	AUC	F1	MCC
86.46 [85.03 – 87.83]	0.806	0.86	0.73	86.46 [85.03 – 87.83]	0.776	0.86	0.73	86.46 [85.03 – 87.83]	0.749	0.86	0.73

Field characteristics within voxels predictive of treatment responders

Within the voxels that discriminate tDCS responders from non-responders (Fig. 4), responders were found to have greater current intensity within these regions (Fig. 5A–B), with greater median (r = 0:811, p < 0:001) and mean current intensity (r = 0:720, p = 0:004) correlated with treatment response (Fig. 5C). Responders produced an effect size of 5.63 (Hedges’ g; Fig. 5D) with a 95% confidence interval between [CI: 3.82–7.94] and a significant two-sided permutation t-test (p < 0:001, 5000 permutations). Behavior change was also related to the azimuthal angle, the angle of the current vector in the axial plane between the electrodes, φ (r = 0:774, p = 0:001) and the x-magnitude of the current vector, J (r = 0:832, p < 0:001). Behavior change was not related to zenith angle, θ (r = 0:289, p = 0:32; see Supplemental Fig. 4), y-magnitude of the current vector, J (r = 0:222, p = 0:45), or the z-magnitude of the current vector, J (r = 0:281, p = 0:33). Thus, current angled toward the cathode related to positive outcomes (see Fig. 6A–C for representative model). On average, the current direction computed within the xy-plane, φ, showed greater percentages of positive angles in responders, whereas non-responders demonstrated greater percentages of negative angles. Thus, an enhanced convergence of current direction toward the cathode was found within responders (Fig. 6D). Shifts in electrode location were also inversely correlated with behavior change (anode: r = − 0:732, p = 0:003, cathode: r = − 0:775, p = 0:001), as well as changes in median current intensity (anode: r = − 0:523, p = 0:06, cathode: r = − 0:632, p = 0:02) and azimuthal angle (anode: r = − 0:579, p = 0:03, cathode: r = − 0:794, p < 0:001) within these voxels.

Fig. 4.

Discrimination maps of regions that predict working memory improvements with the percent contribution of each voxel to the SVM decision function, superimposed onto the MNI152 Template.

Fig. 5.

Plots to demonstrate the current density characteristics within regions predictive of tDCS responders. A) Histogram of current intensity (bin width of 0.0013 Am-2) with the y-axis representing the number of observations in each bin divided by the total number of observations, where the sum of all bar heights is equal to 1. B) Cumulative histogram of current intensity with the height of each step equal to the cumulative number of observations in the bin over the total number of observations in each bin and all previous bins where the height of the last bar is equal to 1. C) Scatter plot of behavior change (post – pre intervention working memory performance) vs. median current intensity. D) The Hedges’ g between responders and non-responders is shown in a Gardner-Altman estimation plot. The mean difference is plotted on a floating axis as a bootstrap sampling distribution. The mean difference is depicted as a dot; the 95% confidence interval is indicated by the ends of the vertical error bar.

Fig. 6.

Single representative participant A) Coronal, B) Sagittal, and C) Axial image of current intensity represented by the color of the images, and current direction represented by the arrows within the images. The color bar represents electric field in volts per meter (v/m). D) Histogram of azimuthal angle δ with the height of each bar representing the number of observations in each bin divided by the total number of observations, where the sum of the all bar heights is equal to 1.

Evaluation of simplified SVM models.

To determine whether the complexity of the original model presented in the current paper is required to achieve sufficient model performance, we assessed a set of simplified SVM models. As electrode displacement was strongly associated with behavior change, we trained an SVM to classify responder and non-responders based on only electrode displacement. In addition, prior research argues that tDCS primarily impacts only gray matter regions. As such, we trained an SVM to classify reponders and non-responders using only current characteristics from gray matter regions. These models did not outperform the original model (see Supplemental Table 1 and Supplemental Table 2; respectively).

Regional contributions toward predictions of tDCS response

Fig. 7A–B illustrates the top ten ranked ROIs based on the mean percent contribution per voxel within each region. See Fig. 7C for the distribution of percent contribution across Harvard-Oxford ROIs. The top ten ROIs that predicted working memory improvements were largely located in the frontal region underneath and between the electrodes (Fig. 7A). These ROIs were labelled as the: 1) Right Superior Frontal Gyrus, 2) Left Superior Frontal Gyrus, 3) Right Middle Frontal Gyrus, 4) Left Putamen, 5) Right Frontal Pole, 6) Right Precentral Gyrus, 7) Left Frontal Pole, 8) Right Pars Opercularis, 9) Left Caudate, 10) Right Supramarginal Gyrus, anterior division (Fig. 7B).

Fig. 7.

A) Visualization of the top 10 regions of interest from the Harvard-Oxford atlas ranked based on their contribution toward predictions of treatment response. B) Rank, label, and mean percent contribution per voxel of the top ten regions of interest. C) a bar graph to represent the average percent contribution per voxel within each ROI of Harvard-Oxford atlas with the top 10 regions of interest highlighted in red.

Dosing charactersitics in responders and non-responders

Equation (6) was used to determine the optimal dosing parameters that are likely to convert non-responders into responders. To match the mean profile of current characteristics demonstrated in responders, non-responders would require an average increase of applied current intensity to 4.3 mA [range(sd): 3.19–5.37 (0.71) mA] and shift in the cathode location 1.25 cm [range(sd): 0.52–2.22 (0.59) cm] closer to its ideal 10–20 location (F3). Fig. 8A and B demonstrate the relationship between differences in electrode placement and behavioral response. Fig. 8C and D demonstrate the required shifts in intensity (8C) and electrode placement (8D) to match the mean profile of responders.

Fig. 8.

A) Anode and B) Cathode displacement between responders and non-responders, with means plotted as bars. Non-responders were found to have greater displacement of both electrodes (anode: F 2 = 6.73, P = 0.023, cathode: F 2 = 19.39, p < 0.001) from their ideal placement compared to responders. Linear regression of C) current intensity and D) current direction, based on the percent difference of average current intensity and cathode displacement versus behavioral change. The optimal stimulation parameters are represented by the diamond in the figure (i.e., the mean within the responder group). * represents p < 0.05 and ** represents p < 0.01).

Discussion

The present study investigated the essential characteristics of tDCS current (i.e., current intensity, current direction, etc.) by using a machine learning algorithm to predict tDCS effects on measured behavioral outcomes. Overall, both current direction and intensity are demonstrated to be critical components of stimulation response. Consistent with our primary hypothesis, the electric field components produced predictions of tDCS response classification beyond chance (i.e., AUC > 50%). While current intensity did marginally outperform other variables, this difference was not statistically significant. Contrary to our second hypothesis, all tested variables were comparable in classifying responders and non-responders. Considering each of these datatypes are derived from the same data vector, it is likely that the interaction of direction and intensity does not provide sufficient new information to outperform intensity alone. Since the SVM relies on pattern similarities between observations to make predictions, the interaction term of direction and intensity did not explain a significant level of new variance that was not already explained by the separate main effects of intensity and direction. Using a small clinical trial dataset, the machine-learning algorithm presented in this paper provides proof of concept that an SVM was able to classify tDCS responders and non-responders with 86% accuracy based on patterns of current characteristics.

Characteristics of current

The weights produced by the SVM algorithm revealed the brain regions that contributed most to predictions of treatment response category. Within these regions, median and mean current intensity were found to positively relate to behavioral response, suggesting greater current intensity may produce greater behavioral response. Since tDCS is typically applied at a fixed dose across individuals, individual differences in anatomy are likely to cause varying amounts of current intensity within these essential brain regions. For instance, those with greater degrees of brain atrophy would have a decreased level of current intensity within stimulated brain tissues and thus may experience reduced efficacy from fixed dose tDCS. It is important to note that while our data demonstrate an association between delivering larger amounts of current intensity and better behavioral responses, this does not necessarily mean that “more is better” universally. For instance, Samani et al. previously reported that applied current doses beyond 2 mA demonstrates nonlinear alterations in neuronal excitability [70]. Thus, within the range of current intensities inducible in the brains of older adults receiving fixed 2 mA dose, increased current intensity appears beneficial for treatment response. Applied current doses beyond this range require further research. Electrode orientation and location during tDCS can also dictate the shape and location of stimulation contact area, altering the pattern of current flow within the brain [22]. Both anode and cathode displacement were found to negatively affect intensity and direction of current within the brain as well as behavioral response. Minor electrode shifts (i.e., ≥1 cm) have been previously demonstrated to dramatically alter the current intensity by up to 38% [21,71]. Therefore, monitoring and increasing electrode placement accuracy via 3D capture techniques [26] may improve individual treatment response.

Regions of importance

The top ten ROIs that predict working memory improvements included the dorsolateral prefrontal cortex (DLPFC), ventrolateral prefrontal cortex (VLPFC), basal ganglia, and cingulo-opercular network regions [72]. The DLPFC and VLPFC are critically involved in monitoring, maintaining, and manipulating task-relevant information. These processes are vital for working memory function [73-75]. Recent studies have reported increased functional connectivity within these regions paired with improvements in working memory performance following applications of tDCS [49,76].

The basal ganglia also play an important role in learning and memory [77]. Specifically, the basal ganglia have been suggested to work conjunctively with the middle frontal gyrus to select information to be stored in working memory [78]. In addition, increased functional connectivity of the cingulo-opercular network (also referred to as the ventral attention/salience network) is associated with better performance on measures of fluid cognition (e.g., executive function) in older adults [79]. The SVM identified critical current features in brain regions previously associated with working memory performance and age-related cognitive decline, serving as an additional proof of principle for this approach.

Future directions

The generated SVM model was highly accurate at classifying treatment responders in older adults, based only on the distributed patterns of electric current. Using these methods, precision tDCS dosing tailored to each individual can be generated to efficiently deliver equivalent current intensity within cortical regions that are essential for producing improvements in working memory. With optimization, application of these presented methods could potentially be expanded to improve tDCS efficiency in not only healthy older adults but also various mental health [7,80,81] and brain-based pathologies [8,9].

Study limitations

The limited number of data points (n = 14) in the current study may affect the generalizability of the model. Thus in our future work, we will use a larger and more heterogeneous tDCS clinical trial dataset (n = 160) [82] (NCT02851511) which is near completion. To maximize the use of available data points and avoid overfitting in this study, we used two-level nested cross-validation to increase the number independent test samples and used 10 permutations to assess the retest reliability of these models. Average accuracy and confidence intervals across these permutations were used to estimate model performance on novel datasets. Our results show a 95% confidence interval of [85.03–87.83%] when the average accuracy is 86.42%, which demonstrates robust performance and small standard deviation of the classification performance. However, cross-validated analyses could lead to overoptimistic results. While efforts have been made to avoid overestimating the performance due to “double-dipping” in the choice of model hyperparameters [83] through nested cross-validation, results should be validated on an independent, larger cohort.

For simplicity of interpretation, we used a binary SVM to distinguish discrete classes: responders and non-responders. However, it is also feasible to utilize machine-learning algorithms to make continuous predictions (i.e., magnitude of improvement). It should also be noted that misclassification of tissue continues to be a challenge for automated segmentation routines [84] which can influence subsequent FEM results since current density estimation is directly related to tissue organization and their assigned conductivity values. For instance, both bone and CSF compartments appear dark in T1-weighted images and thus bone tissue might be misclassified as CSF, or vice versa. This type of tissue misclassification can alter the amount of current entering the skull cavity [85] and thus affecting the estimated current found in the brain that was used for SVM classification. Therefore, the findings of the present study are limited by segmentation inaccuracy produced from the automated segmentation routines. Future studies are warranted to implement manual corrections to automated tissue segmentation and further analyze SVM results for automated versus manual segmentation process. Further, the present study serves as a proof of concept study for classification of tDCS responders and non-responders by utilizing features submitted to the classifier from the computational FEM. Introducing additional neuroimaging modalities and clinical data into the model may further enhance performance and predictive value of machine-learning based approaches for understanding tDCS treatment response.

Conclusion

Clinical applications of tDCS have grown exponentially, yet reproducibility of these studies remains a challenge. At present, the optimal stimulation parameters remain unknown, making it difficult to ensure optimal treatment response on an individual basis. In this study, we proposed novel methods using FEM and SVM algorithms to detect the critical features of tDCS current, and identify stimulated cortical structures implicated in producing intended functional outcomes. We tested our methods in fourteen tDCS participants that underwent multiple days of stimulation at the F3–F4 locations. Our results demonstrated that current intensity has the strongest predictive value for treatment responders, with the performance of current direction only slightly inferior. These data suggest that reducing electrode displacement and delivering greater current intensity to essential brain regions implicated in treatment response are important for producing positive alterations in working memory performance. Future studies may include a larger cohort to generate more generalizable predictive models. Findings from this study can be coupled with customized electrode montages and dosing parameters to potentially enhance functional gains from tDCS treatment on an individual basis for a variety of therapeutic applications.