Dose-response in modulating brain function with transcranial direct current stimulation: From local to network levels

Understanding the dose-response relationship is crucial in studying the effects of brain stimulation techniques, such as transcranial direct current stimulation (tDCS). The dose-response relationship refers to the relationship between the received stimulation dose and the resulting response, which can be described as a function of the dose at various levels, including single/multiple neurons, clusters, regions, or networks. Here, we are focused on the received stimulation dose obtained from computational head models and brain responses which are quantified by functional magnetic resonance imaging (fMRI) data. In this randomized, triple-blind, sham-controlled clinical trial, we recruited sixty participants with methamphetamine use disorders (MUDs) as a sample clinical population who were randomly assigned to receive either sham or active tDCS. Structural and functional MRI data, including high-resolution T1 and T2-weighted MRI, resting-state functional MRI, and a methamphetamine cue-reactivity task fMRI, were acquired before and after tDCS. Individual head models were generated using the T1 and T2-weighted MRI data to simulate electric fields. In a linear approach, we investigated the associations between electric fields (received dose) and changes in brain function (response) at four different levels: voxel level, regional level (using atlas-based parcellation), cluster level (identifying active clusters), and network level (task-based functional connectivity). At the voxel level, regional level, and cluster level, no FDR-corrected significant correlation was observed between changes in functional activity and electric fields. However, at the network level, a significant positive correlation was found between frontoparietal connectivity and the electric field at the frontopolar stimulation site (r = 0.42, p corrected = 0.02; medium effect size). Our proposed pipeline offers a methodological framework for analyzing tDCS effects by exploring dose-response relationships at different levels, enabling a direct link between electric field variability and the neural response to tDCS. The results indicate that network-based analysis provides valuable insights into the dependency of tDCS neuromodulatory effects on the individual’s regional current dose. Integration of dose-response relationships can inform dose optimization, customization, or the extraction of predictive/treatment-response biomarkers in future brain stimulation studies.


Introduction
Non-invasive brain stimulation techniques, particularly transcranial direct current stimulation (tDCS), have gained prominence in investigating the relationship between modulated brain regions and stimulation outcomes [1,2].tDCS has shown promising results in modulating brain activity and connectivity in both healthy individuals and those with neurological or psychiatric disorders [3][4][5].Previous studies have demonstrated that tDCS can alter cortical excitability and brain functions [2,6].However, the inter-individual variability of tDCS poses challenges in detecting intervention effects at the group level [7], and small effect sizes further limit its potential efficacy [8][9][10].
The variability in tDCS responses can be attributed to anatomical and functional factors, including individual differences in fat thickness, skull thickness, cerebrospinal fluid volume, ongoing brain activity, and connectivity [7,[11][12][13].Computational models are commonly employed to estimate the received stimulation dose, as it is not feasible to directly measure cortical electric field (EF) strength invasively [14].These models have been validated through intracranial recordings, physiological measurements, and neuroimaging studies [15][16][17][18][19]; (e.g., see [20,21] for more details on SimNIBS validation and its comparison with other computational modeling method which are commonly used for creating computational head models).Computational head models based on high-resolution structural MRI provide insights into the distribution of EFs across the cortex, with the assumption that EF intensity relates to functional responses at the cortical target site.Investigating the dose-response relationship between the EFs and underlying brain functions can thus help explain the effects of tDCS [22].
To investigate the dose-response relationship and understand how tDCS-induced EFs influence brain functions, multiple modalities of magnetic resonance imaging (MRI) data, including structural and functional MRI, could be informative.Structural MRI is used to create individualized computational head models that predict the distribution patterns of EFs

Summary
In this study, we sought to understand how varying received doses of transcranial direct current stimulation (tDCS), a non-invasive brain stimulation technique, influence brain function as measured by functional magnetic resonance imaging (fMRI).We explored this "dose-response" relationship in participants with methamphetamine use disorders, using them as a representative clinical sample.Leveraging computational models based on individual MRI scans, we predicted cortical electric fields (EFs) as indicators of the received stimulation dose for each participant.Our study aimed to uncover the relationship between these EFs and changes in brain function at various levels: from individual brain locations (voxels) to broader network connections.While no significant correlation emerged between EFs and brain activity when focusing on individual voxels, specific regions of interest, or clusters, a notable correlation was observed at the broader network level, particularly between EFs and the frontoparietal connectivity during drug cue exposure.These results suggest that the effects of tDCS might be more effectively understood by examining wider brain networks rather than just individual points or regions.This research sheds light on the potential of individualized treatment plans using tDCS, underscoring the significance of comprehending its impact on interconnected brain networks.over the cortex.Functional MRI captures cortical functional activity in response to the injected current.Few studies have integrated computational head models with functional MRI data to explore dose-response relationships.
Based on our systematic review, as of the end of Jun 2023, only twelve studies have integrated computational head models with fMRI data to measure dose-response associations in both healthy and clinical populations [1,[22][23][24][25][26][27][28][29][30][31][32].Available dose-response studies could be categorized into two main approaches: regional-based and network-based methods that used a simple linear correlation/regression analysis to quantify dose-response relationships.Previous investigations of dose-response relationships have primarily focused on regional-based association which has been performed at different levels including voxel-level, region-of-interest (ROI)-level, or cluster-level, and some significant (e.g., [1]) or non-significant associations (e.g., [32]) were also reported.However, brain regions are functionally connected, and exploring the relationships between dose (EFs) and response (brain functions measured with fMRI) within brain circuits or networks is crucial.Three published papers have explored the doseresponse at the network-level [26,31,32].For example, it has been shown that current density in the left DLPFC positively correlated with changes in functional connectivity between two predefined ROIs (left dorsolateral and left ventrolateral PFC) during a working memory task, using psychophysiological interaction analysis and simulating precise computational head models [26].Despite the available evidence according to previously published papers, there is no clue about which level of integration from voxel to network would be more informative in dose-response relationship analysis based on integrating head models with fMRI data.
This study aimed to examine the association between individualized received dose (EFs estimated using subject-specific computational head models) and brain response (changes in fMRI data) at multiple levels.Although these four levels of analysis are interconnected, we explored all four levels to determine which approach may be more effective for future doseresponse analyses.By investigating the associations at different levels, we aimed to gain a comprehensive understanding of the dose-response relationship in tDCS and identify the most informative pipeline for studying the effects of stimulation.This analysis allows us to evaluate the relative strengths and limitations of each level, providing insights into which level of analysis may be more effective in elucidating the underlying mechanisms of tDCS-induced brain modulation.Given the highly interconnected nature of the brain, we hypothesize that associations between EFs and functional changes will be stronger at the network level compared to other levels of analysis.Here, we present the data from our pre-registered trial (NCT03382379) [33], which investigated the dose-response relationship in a group of individuals with methamphetamine use disorders (MUDs) as a sample of a clinical population offering novel insights into the effects of tDCS and its potential application in personalized treatment approaches.

Ethics statement
This study was approved by the Western Institutional Review Board (WIRB Protocol #20171742, Transcranial Direct Current Stimulation to Modulate Top-Down Regulation for Drug Craving in Methamphetamine Use Disorder [NeuroMethDC]).The research adhered to the principles of the Declaration of Helsinki and was conducted in compliance with relevant guidelines and regulations.Formal written consent was obtained from all participants before any procedures took place.All participants were adults; thus, no parental or guardian consent was required.

Participants
Participants included 60 individuals (all-male, mean age ± standard deviation (SD) = 35.86± 8.47 years, ranging from 20 to 55) diagnosed with methamphetamine use disorder (MUD).They were recruited during early abstinence from the 12&12 residential drug addiction treatment center in Tulsa, Oklahoma as part of a clinical trial investigating the effectiveness of tDCS in reducing methamphetamine craving (ClinicalTrials.govIdentifier: NCT03382379).Further details regarding the inclusion and exclusion criteria can be found in the supplementary materials S1 Text.

Data collection procedure
The data acquisition procedure for this study is depicted in Fig 1 .It was a randomized, tripleblind, sham-controlled clinical trial with two parallel arms.Participants provided written consent and were randomly assigned to receive either active or sham stimulation.Demographic and substance use profiles of each group and the primary outcomes of the trials are presented in our previously published paper [33].Neuroimaging data, including structural MRI, restingstate, and task-based fMRI, were collected in a pre-stimulation/post-stimulation design.Methamphetamine craving was assessed using a Visual Analogue Scale (VAS) at multiple time points.

Dose-response relationships in tDCS studies
The methamphetamine cue reactivity (MCR) task was conducted following each restingstate scan.The MCR task involved presenting pictorial methamphetamine cues in a block design, consisting of two sets of distinct but equivalent pictures that had been validated in a previous study.The total duration of the task was approximately 6.5 minutes and included 4 blocks of neutral pictures and 4 blocks of methamphetamine pictures.Each block consisted of 6 pictures from the same category (meth or neutral), with each picture presented for 5 seconds and a 0.2-second blank interval between them.A visual fixation point was displayed for 8 to 12 seconds between each block.After each MRI session participants were asked to rate their current level of methamphetamine craving on a scale from 1 to 100 using VAS score.Craving scores were also collected at the baseline session and the day after the stimulation.More details on imaging parameters can be found in supplementary materials S2 Text.

Transcranial direct current stimulation
tDCS was administered using two saline-soaked surface sponge electrodes, each with an area of 5x7 cm2, connected to the battery-driven NeuroConn DC stimulator MR.To target the right DLPFC unilaterally, a bipolar non-balanced configuration was employed based on "Beam F3-system" [34] by taking into account the nasion-to-inion and tragus-to-tragus distances, as well as head circumference, to determine the electrode coordinates for each individual's scalp.The anode electrode was positioned over F4 in the EEG 10-20 standard system, with the long axis of the pad pointing towards the vertex of the head.The cathode electrode was placed over the contralateral eyebrow (Fp1 EEG electrode site, also known as the supraorbital position), with the long axis of the pad parallel to the horizontal plane.To secure the electrodes in place, multiple rubber headbands were used.
During active stimulation, tDCS was delivered for 20 minutes at an intensity of 2 mA, with a 30-second ramp-like fade-in, followed by 19 minutes of active stimulation, and a 30-second ramplike fade-out.For the sham stimulation procedure, the stimulator automatically turned off after 100 seconds, consisting of a 30-second ramp-like fade-in, 40 seconds of active stimulation, and a 30-second ramp-like fade-out.This simulated the typical sensations induced by tDCS.In the sham group, the fade-out was followed by 18.3 minutes without any further stimulation, during which impedance was periodically checked to ensure the average current over time did not exceed 2 μA.Throughout both active and sham tDCS sessions, the impedance was maintained below 10 k-Ohm.

fMRI data preprocessing
Functional data analysis was conducted using the AFNI software package.Preprocessing steps were applied, including the removal of the first three pre-steady state images.The following procedures were performed: despiking, slice timing correction, realignment, transformation to MNI space, and Gaussian FWHM smoothing with a kernel size of 4 mm.To account for potential confounding factors, three polynomial terms and the six motion parameters were regressed out during the analysis.Additionally, TRs with excessive motion were identified and censored based on a criterion defined as the Euclidean norm of the derivative of the six motion parameters exceeding 0.3.
The modeling of the neutral and methamphetamine cue blocks was accomplished by convolving the stimulus timing with a 31-second block regressor, resulting in the estimation of one response for each condition before and after stimulation.Notably, the primary analysis did not incorporate separate regressors for box trials or response periods.

Computational head modeling
High-resolution T1-weighted MR images were used along with SimNIBS 3.1 software to generate personalized computational head models [20].The head models consisted of six main tissues: white matter (WM), gray matter (GM), cerebrospinal fluid (CSF), skull, vitreous bodies of the eyes, and skin.The segmentation process was performed using the "headreco" function in SPM 12 combined with the CAT12 toolbox, and the accuracy of the segmentations was visually inspected against the high-resolution T1-weighted MR images.Tetrahedral volume meshes were created based on the segmented surfaces, with approximately 3 million tetrahedra assigned to each personalized head model.Virtual rectangular pads with dimensions of 5 x 7 and a thickness of 1 mm were modeled and placed on the scalp of each realistic head model, corresponding to the anode and cathode positions over F4 and Fp1, respectively.The conductivity values used for the different tissues in the simulations were based on previously reported values: WM = 0.126 S/m, GM = 0.275 S/m, CSF = 1.654S/m, skull = 0.010 S/m, skin = 0.465 S/ m, and eyeballs = 0.5 S/m [32].To simulate the distribution patterns of the electric field, a current strength of 2 mA was applied, and the electric field (EF = -rφ) was solved using a finite element solver (FEM) under the assumption of a quasi-static regime [35].The tangential electric field, representing the strength of the electric field, and the radial electric field, reflecting the currents entering or leaving the cortex, were calculated for each individual.The simulation results were then transformed into the standard space to enable comparability for group-level analysis.The visualization of the results was performed using Gmsh [36] and MATLAB (version 2019b, The Math Works Inc).General findings about EF simulations can be found in supplementary materials S3 Text.
While SimNIBS is a widely used method for creating computational head models in brain stimulation studies and validated before [20,37,38], it is worth mentioning that there are alternative methods available for indirectly estimating tDCS-induced electric field strength using MRI imaging.One such method is MR-based current density imaging (CDI), which allows for the monitoring of current flow inside the brain during tDCS [39][40][41].CDI involves generating personalized three-dimensional volume conductor models using anatomical MR images and incorporating directional information from diffusion tensor MRI scans (DT-MRI) along with measured magnetic flux density data resulting from the externally applied stimulation current.By assuming an anisotropic conductivity distribution, the method calculates a model-predicted current density and magnetic flux density distribution.Through an iterative process that compares the differences between measured and computed magnetic flux density data, the transversal components of the current density distribution produced by the stimulation current are updated.The CDI technique has been validated using realistic three-dimensional human head models, measured DT-MRI images, and simulated magnetic flux density data, providing reliable and quantitative visualization of internal current density distributions during tDCS [39,40].

Association analysis
A sample size of 30 participants per arm provided 80% power to detect an effect size (Cohen's d) of 0.74.Associations were investigated using Pearson Correlation Coefficients with FDR correction (q was set to 5%).The associations between cortical electric fields (received dose) and changes in brain function (response) were examined at four different linear levels (Fig 2 is a schematic representation of each level).( 1) Voxel-level associations: Analyzing correlations between EFs and fMRI data at the level of individual voxels.Each voxel represents a small three-dimensional element in the brain which includes a complex mixture of neural and nonneural components like neurons, glial cells, blood vessels, and interstitial fluid.(2) ROI-level associations: Assessing changes in brain activation during the drug cue reactivity task by calculating average measures of activity within specific regions of interest (ROIs).ROIs are collections of voxels that are defined based on anatomical criteria such as atlas based parcellation.
(3) Cluster-level association: Investigating the effects of tDCS on whole-brain functional activation by identifying clusters of contiguous voxels that exhibit significant changes in activity.A cluster represents a group of adjacent voxels with similar functional characteristics.(4) Task-based network-level associations: Examining patterns of functional connectivity between a seed region and the rest of the brain or between two seeds during the drug cue reactivity task.Functional connectivity refers to the temporal correlation of activity between different brain regions, and a network represents a set of interconnected regions.Association analyses were performed in the active group.However, due to significant results at the network-level we also checked if there is any significant correlation in the sham group at the network-level.

Regional-based associations 2.8.1. Voxel-level associations.
A block design analysis was conducted on the meth > neutral contrast in pre-and post-stimulation scanning sessions using general linear models (GLM).The functional map and computational head models were transformed to MNI space, excluding non-brain voxels.EFs were resampled to match the resolution of the functional map, resulting in final maps comprising approximately 1050624 (96 x 114 x 96) voxels.Three-dimensional voxels were vectorized, along with the different columns of a matrix for each map, where participants were stacked along the rows.Let Xi and Yi represent the Framework for Data Analysis.The integration of computational head models with fMRI data was achieved through regional-based associations and networkbased associations.Regional-based methods encompassed three levels: voxel-level, regions of interest (ROIs)-level, and cluster-level, while network-based associations involved exploring the relationships between electric fields (EF) and changes in functional activity/connectivity within large-scale networks or functional circuits.https://doi.org/10.1371/journal.pcbi.1011572.g002 column vectors across all participants for the i-th voxel from the EF and fMRI matrices.A whole-brain voxel-wise correlation was calculated between each corresponding voxel within the brain.The spatial distribution of the correlations (r-values and p-values) was then used to reconstruct the brain map.Statistical results were corrected for multiple comparisons using false discovery rate correction (FDR).
2.8.2.ROI-level associations.For ROI-level analysis, brain activation during the drug cue reactivity task was calculated for the meth > neutral contrast before and after stimulation separately.Mean beta weight values were estimated for all extracted ROIs using the Brainnetome atlas parcellation [42] to assess the level of changes in task-related activity during the meth versus neutral condition.Brainnetome atlas was also applied to computational head models in the active groups and averaged EF strength was extracted from each subregion.To determine the relationship between induced EF in each brain region and changes in neural activation (i.e., functional activity after tDCS compared to before stimulation) in response to cue exposure, the correlation between averaged EF strength and changes in functional activity was calculated in each corresponding Brainnetome atlas sub-regions.Then, FDR correction was applied to the statistical results.
2.8.3.Cluster-level association.The effects of tDCS on whole-brain functional activation were studied using a linear mixed-effect model (LME) at the voxel-level, which accounted for the non-independence of observations in repeated measure data.The model included time (pre-and post-stimulation) and group (sham and active), along with their interaction, as fixed effects.The subject was treated as a random effect.The lme function of the R-package nlme was implemented to fit the linear mixed model [43].To identify clusters showing significant changes in functional activation, a cluster-based analysis was performed.This analysis considered both regions of increased activation and regions of decreased activity following tDCS intervention.The group-level functional t-map was computed in MNI space, and family-wise error (FWE) correction was applied using Monte Carlo simulations (3dClustSim, AFNI) with an alpha threshold > 0.1.Results with a significance level of P < 0.005 and cluster size > 40 were considered statistically significant.
For each identified cluster, the averaged BOLD signals were extracted, representing the combined response within the cluster.Cluster masks were then applied to the computational head models to extract the averaged EF intensity from each cluster.Correlations between EFs and the functional map obtained from the LME model were calculated at the cluster level using Pearson Correlation Coefficient analysis.In this study, a cluster-based analysis was employed to offer a comprehensive assessment of the functional changes induced by transcranial direct current stimulation (tDCS).This analytical approach is designed to capture both regions of increased and decreased neural activation, thereby providing a holistic view of tDCS effects on the brain.Importantly, the cluster-based method allows for the detection of distributed effects across various brain regions, encapsulating both positive and negative modulations in functional activity and connectivity.The study specifically focused on temporal changes (pre-vs.post-stimulation) and the interaction between time and experimental groups (active vs. sham) to isolate the effects attributable to the tDCS intervention.While baseline differences between the active and sham groups were possible, these were rigorously accounted for in post hoc analyses to ensure that the observed functional changes were indeed a consequence of the tDCS intervention rather than pre-existing group differences.

Network-based associations
2.9.1.Task-based network-level associations.To examine seed-to-whole-brain taskbased functional connectivity, a seed region was defined based on the location of the maximum EF.The individualized computational head models were transformed into fsaverage space to calculate the averaged EF across the entire population.From the averaged EF map, the 99 th percentile of the maximum EF was determined, and a 10 mm sphere was defined around this threshold.This sphere was used as the seed region for conducting generalized psychophysiological interaction (gPPI) analysis in CONN toolbox [44], which identifies voxels in the brain exhibiting altered connectivity with the seed region during a specific context, such as the drug cue reactivity task in this study.The averaged BOLD signal was extracted from the seed region, and task-modulated functional connectivity was computed.At the first-level analysis, the design matrix included the time course of the seed region, the cue exposure task time course, and the interaction between the task and the BOLD signal in the seed region.For the second-level analysis, an interaction between time (pre-vs.post-stimulation) and group (active vs. sham) was calculated.The voxel-level threshold was set at p uncorrected < 0.01, and the cluster-level threshold was set at p FDR corrected < 0.05.To consider the significant clusters obtained from the seed-to-whole PPI results, the masks of these clusters were utilized.ROI-to-ROI PPI connectivity analyses were then performed between the seed region located around the 99 th percentile and the cluster masks.The correlation between ROI-to-ROI PPI connectivity and averaged EFs in the seed region was calculated in each active and sham group separately.

Considering behavioral data
Behavioral data, specifically VAS scores obtained through the question 'How much craving do you have right now?', were collected immediately before and after each MRI session.Significant associations were observed between EFs and BOLD signal at the network-level, indicating a connection between the measured neural activity and the distribution patterns of EFs.To further investigate the relationship between the obtained results and behavioral outcomes in both the active and sham groups, additional correlations were computed between VAS scores, BOLD signal, and EFs.

Regional-based associations
3.1.1.Voxel-level results.At the whole-brain voxel-level, block-designed analyses revealed no FDR-corrected significant correlation between electric fields and BOLD signal changes (p corrected > 0.05).Only a small proportion of voxels exhibited a small effect size (0.1 < |r| < 0.3) (9.74% of the voxels), a medium effect size (0.3 < |r| < 0.5) (1.36% of the voxels), and a mere 0.09% of voxels showed a large effect size (0.5 < |r| < 1).The total number of voxels analyzed was 96x114x96, and the Pearson Correlation Coefficient (|r|) was used as a measure of effect size (Fig 3).

ROI-level results.
At the ROI-level, pairwise correlations between EFs and changes in fMRI BOLD signal (Post-Pre) was conducted in the meth vs. neutral contrast, using atlasbased parcellation of individualized head models and fMRI maps with the Brainnetome atlas (Fig 4).At the whole-brain ROI-level, results revealed no FDR-corrected significant correlation between EFs and BOLD changes considering p corrected > 0.05.However, without FDR correction, we identified 11 ROIs that showed a correlation with p uncorrected < 0.05.Among these ROIs, a small proportion exhibited a medium effect size (0.3 < |r| < 0.5), accounting for 4.29% of the total ROIs analyzed, while only 0.95% of the ROIs displayed a large effect size (0.5 < |r| < 1) including two brain regions Brodmann area 11 medial area located in orbital gyrus (r = 0.58, p uncorrected = 0.001) and Brodmann area 21 rostral area in middle temporal gyrus (r = -0.55,p uncorrected = 0.002).Only the first region received EFs above 0.2 V/m, however, did not survive FDR correction.

Network-based association
3.2.1.Network-level results.Our gPPI results (Fig 6) showed two significant clusters as follows: (1) located in bilateral amygdala and hippocamp (cluster size: 1356, (x, y, z) peak coordinate in MNI space: (18,-2,-30) and t-value: 6.59), and (2) located in right inferior parietal lobule (IPL) (cluster size: 701, (x, y, z) peak coordinate in MNI space: (46,-38,38), and t-value: 5.20).ROI-to-ROI gPPI connectivity was also calculated between seed region and each cluster and correlation between averaged EF in seed region and seed-to-cluster gPPI connectivity was calculated.For the network-level analysis, frontoparietal connectivity showed a positive significant correlation with EF in the frontal stimulation site (r = 0.43 (medium effect size), p corrected = 0.03).We also checked the correlation results between simulated EFs and frontoparietal task-based connectivity in the sham group and no significant correlation was found (r = 0.10, p corrected = 0.61).

Behavioral results and correlation with EF and neural response at the network-level
Participants rated their cravings based on the Visual Analog Scale (VAS) ranging from 0 to 100 in both active and sham stimulation groups.The linear mixed-effects (LME) analysis revealed a significant effect of time (P < 0.0001), indicating changes in craving scores over time.However, there was no statistically significant effect of group or time × group interaction (Fig A in S4 Text).
During the pre-stimulation phase, both the active and sham groups showed a significant increase in craving scores (P uncorrected < 0.0001 for both groups) after the fMRI cueinduced craving task in the first scan (T2 time point in  assessment (T1).The mean craving scores ± SE were 59.07 ± 6.17 for the sham group and 64.40 ± 5.20 for the active group at T2, while at T1, the scores were 39.4 ± 6.07 for the sham group and 42.50 ± 5.26 for the active group.Our results demonstrated that there was no significant increase (P uncorrected > 0.05) in self-reported craving after the second exposure to drug cues (T4) compared to before the second cue exposure (T3).The mean craving scores ± SE were 24.67 ± 4.96 for the sham group and 27.23 ± 4.63 for the active group at T4, while at T3, the scores were 17.83 ± 21.55 for the sham group and 25.93 ± 3.77 for the active group.Furthermore, a significant reduction in craving scores (P uncorrected < 0.05 for both groups) was observed after stimulation (T3) compared to the before stimulation (T2) for both the sham and active groups.The mean craving scores ± SE were 17.83 ± 21.55 for the sham group and 25.93 ± 3.77 for the active group at T3, while at T2, the scores were 59.07 ± 6.17 for the sham group and 64.40 ± 5.20 for the active group.Additionally, a significant difference in craving scores was found between the baseline (T0) and the day after stimulation (T0: sham 49.37 ± 5.74, active 61.37 ± 4.50; T5: sham 11.17 ± 3.60, active 12.67 ± 3.33).
Despite the absence of significant differences between the active and sham groups in terms of craving self-reports, our exploration at the network level revealed significant correlations between electric fields (EFs) and neural responses only in the active group.To further investigate potential relationships between EFs, neural responses, and behavioral outcomes, we examined the correlation between behavioral data (changes in VAS scores) and neural responses or EFs at the network level.However, our analysis did not identify any significant correlations between changes in behavioral data and either frontoparietal PPI connectivity (meth vs. neutral) changes or EFs at the network level.We only found a significant correlation (r = 0.37, p uncorrected = 0.049) between frontoparietal PPI connectivity during meth cue exposure in the pre-stimulation scan (at first scan which can be named baseline) and changes in craving score (VAS changes from T3 to T4 which can be named second cue-induced craving) after active stimulation.

Discussion
In this study, we explored the impact of electric fields (EFs) on brain function, assessed through fMRI, across four analytical levels: voxel, ROI, cluster, and network.Voxel and ROIlevel analyses showed no significant correlation between EFs and BOLD signal changes after multiple comparisons correction.Cluster-level analysis identified five active clusters in the left hemisphere with a time-by-group interaction effect in response to tDCS, but these did not correlate with EFs.Network-level analysis, however, revealed significant correlation between the frontal stimulation site and frontoparietal networks, specific to the active stimulation group.Despite these neural changes, no significant differences in self-reported craving scores were observed between the active and sham groups.Finally, we found a significant relationship between baseline frontoparietal connectivity and craving score changes post-stimulation, but this did not correlate with EF-induced network changes, suggesting that EF effects on neural activity may not directly influence subjective craving levels.Taken together, our findings highlight the complexity of EF effects on neural and behavioral outcomes, underscoring the need for multi-level analyses to fully understand the impact of electrical stimulation interventions.

Voxel-level integration
The voxel-level analysis in our study did not reveal any significant correlation between EFs and changes in BOLD signal, which aligns with our initial hypothesis.In contrast, the only previous study that investigated the correlation between EFs and fMRI at the voxel-level reported a significant association between EF distribution patterns in a standard brain and a group-level functional map using a voxel-level approach [27].It is important to note that the methodology employed in that previous study differed from ours.They calculated a single correlation across all voxels using a standard head model in MNI space and a group-level functional map derived from second-level analysis of fMRI data [27].In our study, we considered individualized EF and BOLD signals for each participant separately and then calculated the correlation for each single voxel across the population in the active group.These methodological distinctions may contribute to the observed discrepancy in findings.Furthermore, it is worth mentioning that our results focused on BOLD signal changes, while the previous study reported the correlation using cerebral blood flow (CBF) data [27].
The interpretation of voxel-level results in the context of transcranial electrical stimulation poses several challenges.One important consideration is the spatial alignment of the EF and BOLD signal maps.For meaningful correlation or regression analyses, it is crucial to ensure that both maps are in the same space, either subject space or an averaged standard space, with the same resolution [30].This alignment enables the examination of associations between corresponding voxels across the population.Furthermore, voxel-wise integration faces limitations due to the inherently low signal-to-noise ratio (SNR) of individual voxel data, which can impact the reliability of associations [45].The large number of voxels involved in voxel-wise analysis also necessitates rigorous correction for multiple comparisons to avoid spurious findings.Consequently, the statistical power to detect significant associations at this level is generally low.
Despite the lack of significant correlations at the voxel level, it is important to consider the potential localized effects of EFs on neural activity.Although the overall voxel-level analysis did not reveal significant associations, a small proportion of voxels exhibited large effect sizes, suggesting that localized regions may respond more prominently to EF-induced effects.This localized influence may be attributed to the intricate interplay between stimulation parameters, neural circuitry, and individual differences in brain structure and function.Understanding the factors contributing to these localized effects can provide valuable insights into the optimization of transcranial electrical stimulation protocols and individualized treatment approaches.In summary, our voxel-level analysis did not identify significant correlations between EFs and BOLD signal changes across the entire cortex.However, the challenges inherent to voxel-wise analysis, including spatial alignment, SNR limitations, and the need for rigorous multiple comparison correction, should be considered when interpreting these findings.The identification of a small proportion of voxels with large effect sizes suggests potential localized effects of EFs on neural activity.Future studies should further investigate these localized effects and explore the factors influencing the variability in EF-induced responses across different brain regions.

ROI-level
At the ROI-level, we applied Brainnetome atlas parcellation to both individualized head models and functional maps in order to extract averaged EF intensity and beta values.While we observed strong correlations in specific brain regions, such as the orbital gyrus and middle temporal gyrus, none of these regions survived the false discovery rate (FDR) correction, as we considered all brain subregions in the Brainnetome atlas for our analyses.In this regard, previous studies have reported significant correlations within predefined regions of interest; however, there is a lack of consensus regarding rules or hypotheses for ROI selection.For instance, among studies that examined dose-response relationships at the ROI-level, some studies defined spherical seed regions with varying sizes, including 2.5 mm, 6 mm, 8 mm around specific coordinates such as the peak activation voxel in fMRI analysis results, mean center coordinates obtained from MRS analyses, or around the hand knob region [1,22,24].The remaining studies utilized different types of atlas-based parcellations [23,46].Despite the methodological differences, the ability to detect significant correlations at the ROI-level is highly dependent on how ROIs are defined.These definitions can vary widely, from using standard anatomical or functional atlases for parcellation to placing small spheres around specific brain coordinates.The chosen method for defining ROIs plays a critical role in successfully identifying meaningful relationships between dose and neural response.For example, Esmaielpour et al. employed two different approaches for ROI definition: anatomical atlas parcellation of the prefrontal cortex and parcellation based on EF strength [46].Similar to our study, they extracted averaged EF strength and beta values in the general linear model (GLM) analysis of task-based fMRI data.They reported a single significant correlation in the frontopolar cortex, which exhibited the highest EFs among ten brain regions.However, no significant correlation was found based on anatomical atlas parcellation.This study underscored the importance of the parcellation method and ROI definition in dose-response relationship studies [46].
In dose-response relationships using regional-based approaches (including both ROIs or clusters), it is essential to extract features from each region.This involves averaging the estimated parameters, such as beta activation maps and EF intensities, across the voxels within the predefined ROIs to calculate the relationships between fMRI-derived parameters and transcranial electrical stimulation-induced EFs.The ROI-based approach overcomes the low SNR observed in voxel-wise analysis by aggregating data from multiple voxels.Moreover, performing statistical testing on ROIs rather than voxels reduces the number of tests conducted.However, this approach comes with certain limitations.Results in brain regions outside of the predefined and hypothesized areas may be overlooked, and the selection of appropriate brain regions can be challenging, requiring specific knowledge about the research problem for defining relevant areas.
Taken together, our ROI-level analysis demonstrated the absence of significant correlations between EFs and BOLD signal changes when considering multiple comparison corrections.This suggests that the effects of EFs may be region-specific and not uniformly distributed throughout the brain.The main challenge in finding significant correlations at the ROI-level corresponds to the absence of a hypothesis to appropriately define ROIs.Further investigations are needed to better understand the underlying mechanisms and functional implications of these region-specific effects.Additionally, future studies should carefully consider the methodology for defining ROIs and the potential trade-off between hypothesis-driven analysis and the risk of missing significant associations in unexpected brain regions.

Cluster-level
We did not find any significant correlation between EFs and changes in functional activity within the active clusters, despite the significant effect of group observed in all five clusters.Consistent with our findings, there is no previous study that reported cluster-based associations based on integrating head models with fMRI data.As we defined in Fig 2, all voxel-level, ROI-level, and cluster-level correlation is a kind of regional-based association.However, the cluster-based integration approach bridges the gap between voxel-level and ROI-level analyses by considering each voxel's priority for integration while performing the analysis within contiguous clusters of voxels, similar to the ROI-based approach.However, several factors may contribute to cluster-level findings.First, it is possible that the complex relationship between EFs and neural responses within the identified clusters is not adequately captured by the cluster-level linear associates.The cluster-based integration approach relies on the aggregation of voxel-level data within a cluster, potentially averaging out subtle but significant effects present at the voxel level.Furthermore, the sensitivity and specificity of the cluster-level analysis may be influenced by factors such as the size and distribution of the clusters, as well as the statistical thresholds applied.

Network-level
In line with our hypothesis, significant dose-response relationships were observed at the network level, highlighting the importance of examining global network effects rather than local associations found in voxel-level, ROI-level, or cluster-level approaches.Specifically, we observed a significant positive correlation between task-based frontoparietal functional connectivity and EF intensity within the prefrontal cortex.These findings align with previous research conducted by Indahlastria et al. during a working memory task, where they reported a significant correlation between current density values in DLPFC and changes in task-based functional connectivity within the working memory network in older adults [26].Another recent study also examined associations at the network-level and reported that the individualized electric field in cortical regions with regional EFs exceeding 99 th and 98.5 th percentile significantly predicts tES-induced outcomes in terms of functional connectivity within sensorimotor areas in the ipsilateral hemisphere in motor cortex stimulation of twenty-five chronic stroke individuals [31].In another recent study [32], two different approaches were used for extracting EFs, to determine the highest EFs the mask included 1,000 gray matter voxels for which the highest EF values were estimated.Alternatively, to determine EFs in cortical regions that are related to the EEG alpha band activity, the mask included 1,000 occipital gray matter voxels with the strongest negative correlation between alpha amplitude and BOLD signal.They reported a significant correlation between EF in alpha-BOLD correlated region and changes in alpha amplitude, however, no significant correlation was observed between EFs and changes in functional connectivity within the default mode network during resting-state fMRI data among twenty-two healthy subjects [32].
Two main points should be highlighted based on our results.Firstly, our findings demonstrated significant associations exclusively at the network level, indicating that the neurobiological effects of tDCS may not manifest in a straightforward manner at the voxel, region, or cluster level.This emphasizes the significance of examining network-level changes during stimulation.In recent years, there has been a growing emphasis on network-based effects of brain stimulation technologies, and our results further support the notion that tDCS effects are likely to spread across networks rather than being confined to specific local brain regions.Our networklevel analysis sheds light on the significant modulatory effects of EFs on functional connectivity within the brain.These findings support the growing emphasis on network-level changes during tDCS and provide insights into the complexity of EF-induced effects on network dynamics.Further investigations should explore the underlying mechanisms of these network-level effects and consider the broader implications for understanding the functional consequences of tDCS.

Considering SUDs as a clinical population
In this study, our specific focus was on individuals with methamphetamine use disorders (MUDs) as a representative sample of the clinical population.There are some other findings in terms of dose-response relationships in SUDs.For example, a pre-post tDCS-fMRI study with F4-F3 electrode montage (bilateral frontal with anode on the right) conducted on MUD participants revealed a significant correlation between changes in brain functions in response to drug cues and the averaged EF strength within the frontopolar area, as opposed to other prefrontal regions [46].Moreover, our previous investigations have examined the impact of substance use disorders (SUDs) on electric field distribution patterns in relation to anatomical alterations [47].Specifically, we focused on individuals with cannabis use disorders, comparing them to healthy controls within the context of tES studies, and reported significant structural alterations in participants with cannabis use disorders compared to healthy controls [47].Additionally, studies involving transcranial magnetic stimulation (TMS) have explored the influence of such alterations and revealed lower gray matter volume in multiple areas of interest for TMS treatment associated with alcohol use disorders (AUDs) [48].These findings suggest that alterations in brain structure due to SUDs may affect EF distribution patterns and stimulation outcomes as previous studies have shown that tES/TMS-evoked in cortical and subcortical activity depend on brain structure.Further investigations are needed to assess the generalizability of the results obtained in this study to healthy subjects [49].

General challenges in integrating head models with fMRI
Integrating EFs with fMRI data presents several challenges that impact result interpretation and accuracy.First, combining a static map (head model) with a dynamic map (fMRI) is difficult due to the complex interactions between targeted and non-targeted brain regions.This complexity makes it challenging to establish a direct relationship between EFs and changes in brain activity measured by BOLD signals.Interpreting the results of brain integration studies can be complex in part due to the intricate interconnections between brain regions.When analyzing the relationship between EFs and BOLD signals, it's crucial to consider both positive and negative correlations.For instance, higher EF strengths might reduce functional connectivity, while lower strengths could enhance it.This complexity calls into question the specific mechanisms through which EFs influence brain networks.Additionally, it's important to recognize that the commonly used statistical measures of connectivity lack a robust biophysical basis.While there are plausible reasons why higher EFs could reduce connectivity, such as the orientation of EFs and neurons, the concept of connectivity itself is largely statistical.Therefore, further research is needed to elucidate the biophysical foundations of connectivity and its interplay with EFs for a more nuanced understanding of brain stimulation outcomes.
Another challenge is selecting appropriate measures for EF and fMRI data, as different metrics provide distinct information about EF distribution patterns and brain activity.In addition to EF strength, the orientation of the EF plays a role in adjusting the "effective" intensity when it is not aligned with the axis from the dendrites to the axon of a neuron.It can also indirectly indicate which neurons are recruited, as those that are most aligned with the vector [50].However, these considerations are all encompassed within the measurement of EF intensity.The accuracy of the modeling process for computational head models is crucial, as inaccurate modeling can affect the distribution of EFs.Precise tissue segmentation and conductivity assignment are essential for accurate computational head models [51][52][53].The choice of integration space, whether subject space or standard space, also affects the results.While standard space integration allows for comparative analysis, using personalized models based on individual scans yields higher accuracy.Overcoming these challenges is vital to obtain reliable insights into the relationship between EFs and neural activity measured by fMRI.Future research should address these challenges and improve the integration process.

Limitations and future directions
While the current dose-response relationship studies exploring the associations between EFs and brain functions show promise, there are several areas that require further investigation.Integrating head models with fMRI is still a relatively new approach in tES studies, and although previous works have reported novel findings, replication studies are necessary to establish a comprehensive understanding of the relationship between EF distribution patterns and changes in BOLD signals.
The complexity of the underlying biophysical processes and the diverse nature of neural responses make it difficult to directly link electric field distribution patterns to the efficacy of tDCS.The predictive value of integrating EF distribution with fMRI data in determining treatment outcomes requires further investigation.To drive the field forward, several important steps should be considered.First, future research should focus on elucidating the biophysical mechanisms underlying the effects of EF distribution on neural responses [54][55][56].This could involve combining computational modeling, animal studies, and advanced neuroimaging techniques to establish a more comprehensive understanding.Second, efforts should be made to develop standardized protocols and methodologies for integrating EF distribution and fMRI data as a measure of dose-response relationship [57].This includes refining the methods for assessing and quantifying EF distribution, improving the sensitivity and specificity of fMRI measurements, and establishing consensus guidelines for data analysis and interpretation.Third, in this study, the sample size is relatively small.Additional investigation involving a larger sample size and more diverse participants (e.g., including females to study the potential role of sex differences and other types of substance use) is needed to generalize the results of this study.Larger-scale studies involving diverse populations and clinical conditions are needed to enhance the generalizability of the findings.By expanding the scope of the investigation, we can better assess the efficacy of tDCS interventions and identify potential predictors of treatment response.
Theoretically, the polarization of neuronal membranes correlates linearly with the electric field they're exposed to [58][59][60].However, this doesn't mean that neurophysiological, behavioral, or clinical outcomes also follow a linear progression.Advancements in computational resources and machine learning techniques, such as deep learning, offer opportunities to explore more complex nonlinear and higher-order relationships between EFs and functional activity/connectivity at the subject or group level.For example, deep learning models can be trained using computational head models and fMRI data for classification or prediction purposes.Machine learning methods can also help identify responders and non-responders based on EFs in targeted and non-targeted brain regions, as well as the initial brain state.This approach may facilitate the identification of multi-modal biomarkers by leveraging the predictive role of baseline brain activation and individualized EFs.
The emergence of closed-loop tES-fMRI studies, where fMRI and tES are performed simultaneously, opens up exciting possibilities for real-time integration of personalized ongoing brain function with EFs [61].This approach enables the optimization of stimulation parameters based on the interaction between EFs and neuronal activity.However, implementing closed-loop head model-fMRI integration in real-time poses challenges, as the online processing of computational head models and fMRI data is time-consuming [62].Future studies will require high-speed algorithms and powerful computational tools to enable efficient implementation of closed-loop integration.In addition, future studies can consider integrating other modalities such as EEG, DTI, or behavioral data with computational head models and fMRI.Integrating multiple sources of information allows for the exploration of hybrid activation patterns, leading to the identification of putative trait-state relationships that could offer more informative insights than relying on a single data source.

Conclusion
In conclusion, this study investigated the integration of EFs with fMRI data and explored the associated challenges and implications.Overall, our findings underscore the complexity of the relationship between EFs, neural activity, and functional connectivity.The region-specific effects of EFs and their modulation of network-level connectivity provide valuable insights into the targeted and context-dependent nature of tDCS effects.Although we only found FDR-corrected significant results at the network level, we also observed strong correlations at the ROI level.However, the lack of a specific hypothesis for selecting brain regions in the ROIlevel analysis may have impacted its efficacy, while the network-based approach was more contingent as it helped to limit our search space and identify significant results.
Nevertheless, the absence of significant correlations between EFs and behavioral outcomes highlights the need for further investigation to fully understand how EF-induced neural changes translate into subjective experiences and behaviors.Future research should delve deeper into the underlying mechanisms of these relationships, explore additional analytical approaches, and consider other factors that may influence the effects of EFs on brain function.By addressing these avenues, we can advance our knowledge and maximize the potential of EF-based interventions for both clinical and cognitive applications.

Fig 1 .
Fig 1.Data acquisition procedure.60 participants with methamphetamine use disorders were randomly assigned to the active or sham group.MRI data including structural MRI for creating head models, resting-state fMRI, and fMRI drug cue reactivity task were collected before and after unilateral DLPFC stimulation.Behavioral data including VAS scores were also collected immediately before and after MRI sessions, at the baseline, and the day after.Abbreviations: sMRI: structural MRI, rs-fMRI: resting-state fMRI, ts-fMRI: task-based fMRI, EF: electric field.https://doi.org/10.1371/journal.pcbi.1011572.g001

Fig 2 .
Fig 2.Framework for Data Analysis.The integration of computational head models with fMRI data was achieved through regional-based associations and networkbased associations.Regional-based methods encompassed three levels: voxel-level, regions of interest (ROIs)-level, and cluster-level, while network-based associations involved exploring the relationships between electric fields (EF) and changes in functional activity/connectivity within large-scale networks or functional circuits.

Fig 3 .
Fig 3. Voxel-level associations.I. Whole-brain voxel-level correlation analysis pipeline.Dose-Response (EF-fMRI correlation) at the voxel-level in a wholebrain approach.a. Electric field (EF) maps were vectorized and arranged along the columns, with participants in the active group represented along the rows to construct the EF matrix.A similar approach was employed to generate an fMRI matrix for the meth > neutral contrast.The active group comprised 30 participants, and the total number of voxels analyzed was 96 x 114 x 96.The correlation vector captures the relationship between each corresponding voxel in the EF and fMRI matrices.II.Whole-brain voxel-level correlation analysis results.a. Group-level functional map.Changes in BOLD signal (Post-Pre) in the meth vs. neutral contrast at the group-level for participants in the active group.b.Group-level analysis of head models.EF distribution patterns observed in participants in the active group, with anode/cathode positioned over F4/Fp1.c. Correlation results.Prior to the false discovery rate (FDR) correction, some voxels showed |r|>0.3 and p uncorrected< 0.05.However, none of them survived FDR correction.https://doi.org/10.1371/journal.pcbi.1011572.g003

Fig 1 )
compared to the initial

Fig 4 .Fig 5 .Fig 6 .
Fig 4. ROI-level associations.a. Cortical BOLD signal change.Applying Brainnetome atlas parcellation to the first-level functional maps obtained from general linear models (Post minus pre-stimulation and meth vs. neutral contrast) for the active group.b.Individualized computational head models.Head models were parcellated using the Brainnetome atlas.Averaged electric fields (EFs) were calculated only in cortical subregions.c.Correlation r values.Pearson correlation coefficients between all pairs of ROIs (EFs and BOLD signal change) were calculated in the active group and none of them survived FDR correction.Bars show mean value and error bars show standard errors.https://doi.org/10.1371/journal.pcbi.1011572.g004