Group and individual level variations between symmetric and asymmetric DLPFC montages for tDCS over large scale brain network nodes

Two challenges to optimizing transcranial direct current stimulation (tDCS) are selecting between, often similar, electrode montages and accounting for inter-individual differences in response. These two factors are related by how tDCS montage determines current flow through the brain considered across or within individuals. MRI-based computational head models (CHMs) predict how brain anatomy determines electric field (EF) patterns for a given tDCS montage. Because conventional tDCS produces diffuse brain current flow, stimulation outcomes may be understood as modulation of global networks. Therefore, we developed a network-led, rather than region-led, approach. We specifically considered two common “frontal” tDCS montages that nominally target the dorsolateral prefrontal cortex; asymmetric “unilateral” (anode/cathode: F4/Fp1) and symmetric “bilateral” (F4/F3) electrode montages. CHMs of 66 participants were constructed. We showed that cathode location significantly affects EFs in the limbic network. Furthermore, using a finer parcellation of large-scale networks, we found significant differences in some of the main nodes within a network, even if there is no difference at the network level. This study generally demonstrates a methodology for considering the components of large-scale networks in CHMs instead of targeting a single region and specifically provides insight into how symmetric vs asymmetric frontal tDCS may differentially modulate networks across a population.

www.nature.com/scientificreports/ any given brain regions as currently pursued in modeling studies, the network-led approach for EF comparison in the brain may represent a more productive strategy when comparing stimulation montages. In a recent study 37 ,Fischer et al. suggested that the results of stimulating a brain target may be strengthened or weakened by other brain regions that work together with the targeted region as a network. They hypothesized that modulation of one brain region may impact and be impacted by other regions through the networks and designed a novel multifocal electrode array at the subject-level that stimulate left M1 with excitatory effects, but simultaneously inhibit activity in other remaining nodes of the motor cortex network. Their network approach for placement of the electrodes showed a greater impact on M1 excitability than stimulating just left M1 alone 37 . In another study, through a group-level analysis of computational models, Gomez and his colleagues used a functional atlas for parcellation of the cerebellum. They investigated how the current spreads in the cerebellum networks and reported significant differences between various electrode arrangements for stimulating the cerebellum 38 . However, none of the studies have taken into account the anatomical differences among subjects' brain variability yet and it remains unclear how much the electrode location can change EF distribution patterns across nodes of large-scale functional networks in a group of subjects with different anatomy.
To bridge this gap, in the present study, we advance a systematic pipeline to investigate EFs inside the brain networks in a group of subjects. We have used structural data from a population with methamphetamine use disorders (MUDs) for generating CHMs. Two of the most common electrode montages for modulating the DLPFC were simulated 39 ; an asymmetric "unilateral" montage and a symmetric "bilateral" montage. Our approach allows identifying affected nodes in large-scale functional networks in the standard space and will be insightful for developing hypotheses on the network modulatory effects we should expect. The results of this study provide a better understanding of the network-based EF distribution patterns of the montages targeting DLPFC and can be used for guiding the selection of electrode arrangements for future network-based modulation of the brain.

Methods and materials
Participants and data acquisition. Participants included 66 subjects (all-male, mean age ± standard deviation (SD) = 35.86 ± 8.47 years ranges from 20 to 55) with methamphetamine use disorder (MUD). All subjects were recruited during their early abstinence from the 12&12 residential drug addiction treatment center in Tulsa, Oklahoma in the process of a clinical trial to measure the efficacy of tDCS on methamphetamine craving (ClinicalTrials.gov Identifier: NCT03382379). None of the participants had a history of neurological disorders, head injury, or other abnormalities demonstrated by sMRI. Written informed consent was obtained from all participants prior to the scans and the study was approved by the Western IRB (WIRB Protocol #20171742). After receiving each subject's written consent, structural MRIs were obtained on a GE MRI 750 3T scanner. High resolution structural images were acquired through magnetization-prepared rapid acquisition with gradientecho (MPRAGE) sequence using the following parameters: TR/TE = 5/2.012 ms, FOV/slice = 24 × 192/0.9 mm, 256 × 256 matrix producing 0.938 × 0.9 mm voxels and 186 axial slices for T1-weighted images and TR/ TE = 8108/137.728 ms, FOV/slice = 240/2 mm, 512 × 512 matrix producing 0.469 × 0.469 × 2 mm voxels and 80 coronal slices for T2-weighted images. Gyri-precise CHMs were generated from a combination of high-resolution T1-and T2-weighted MR images for 52 subjects. For a subset of participants (N = 14), T2-weighted MRIs were not available and head models were created only based on T1 images. This study was conducted in accordance with the Declaration of Helsinki and all methods were carried out in accordance with relevant guidelines and regulations.
Creation of head models and tDCS simulation. The flow diagram of the modeling and data extraction pipeline is briefly shown in Fig. 1. Head models were generated for all of the subjects to visualize how current flows through the brain in each individual. Modeling of the EF distribution patterns was performed using the standard SimNIBS 3.08 pipeline 40 . Automated tissue segmentation was performed in SPM 12 combined with CAT12 toolbox via the SimNIBS "headreco" function. The head volume was assigned to six major head tissues including white matter (WM), gray matter (GM), cerebrospinal fluid (CSF), skull, scalp, and eyeballs. Segmentation results were visually examined slice-by-slice to ensure correct tissue classifications. Segmented surfaces were used to create tetrahedral volume meshes and about 3 × 10 6 elements were assigned to each head model.
Two of the most widely studied tDCS conventional electrode montages for targeting DLPFC were simulated; bipolar asymmetric and bipolar symmetric arrangements 41 . For both montages, 5 cm × 7 cm rectangular pads with 1 mm thickness were generated virtually. In the first montage, as a bipolar asymmetric arrangement (F4-Fp1), the anode electrode was positioned over the F4 electrode location on the standard 10-20 EEG international system with the long axis of the pad pointing towards the vertex of the head. The cathode electrode was positioned over the contralateral eyebrow (Fp1 EEG electrode site also referred to as the supraorbital position) with the long axis of the pad parallel to the horizontal plane. In the second montage, symmetric bipolar (F4-F3), anode and cathode electrodes were respectively placed over F4 and F3 electrode locations on EEG standard system with the long axis of the pads pointing towards the vertex. The locations and directions of the anodes are thus similar in both configurations and the only difference between montages is the position of the cathode electrodes. In order to have the same electrode placement for all participants, we manually placed electrodes in the SimNIBS simulation window to precisely check the electrode location and orientation for each individual. For the electrode location, we used the standard EEG cap (EEG10-10_UI_Jurak_2007) for each subject to place electrodes over F3, F4, or Fp1. In order to have a similar electrode orientation among the population, we visually checked the orientation for each individual. We visualized simulation results in "Gmsh" software, one by one, to make sure that all subjects have similar electrode placement in terms of electrode location and orientation. These arrangements are also referred to as unilateral and bilateral montages in DLPFC stimulation 42,43 under the assumption the cathode is functionally inert for the asymmetric (unilateral) case. www.nature.com/scientificreports/ After electrode placement, the simulations were run with a current amplitude of 2 mA and EFs were calculated based on the finite element method (FEM). To generate CHMs, electrical conductivities were assumed to be constant. The assigned isotropic conductivities were based on the previously reported values; WM = 0.126, GM = 0.275, CSF = 1.654, skull = 0.010, skin = 0.465, eyeballs = 0.5 all in S/m (siemens/meter) 44 . The results were visualized using Gmsh 45 and MATLAB (version 2019b, The Math Works, Inc.). The absolute values of the EFs (|EF|) were calculated at each brain node. In order to make the individualized EF maps comparable, the simulation results for all subjects/montages were transformed from the original native space to the standard average surface ('fsaverage') of FreeSurfer (http://surfe r.nmr.mgh.harva rd.edu). In this way, the same coordinates in different subjects correspond to the same location in CHMs. This normalization makes statistics possible at the group-level and allowing for the comparison of the EF distribution patterns between two montages. Data analysis. Functional atlas parcellation was used for comparing EFs in two montages as distributed over large-scale brain networks. Yeo7-2011 surface-based atlas was used for extracting the topology of the seven resting-state large-scale networks including visual (Vis), somatomotor system (SomMot), dorsal attention (Dor-sAttn), ventral attention (VentAttn), limbic system (Limbic), executive control network (ECN), and default mode network (DMN) 46 . Since each parcellated network in Yeo7 atlas covered a widely distributed area, we applied Schaefer-100 atlas to CHMs for extracting the finer sub-regions of each functional network 47 to determine the amount of current reaching different parts within a network. We then merged the areas that were placed adjacent to each other to form the main nodes of the networks. Left VentAttn was divided into four main nodes including tempro-occipital-parietal (TempOccPar), frontal-operculum-insula (FrOperIns), prefrontal cortex (PFC; DLPFC node), and medial (Med) nodes. Right VentAttn network parcellated into three sub-networks including TempOccPar, FrOperIns, and Med nodes. The limbic network has two main nodes in each hemisphere including orbito-frontal cortex (Limbic-OFC) and temporal pole (Limbic-TempPole) nodes. The main nodes in the left and right ECN were tempro-parietal (ECN-TempPar), ECN-PFC, and precuneus cingulate (ECN-PcunCing). DMN was also divided into three main nodes in each hemisphere including DMN-TempPar, DMN-PFC, and precuneus posterior cingulate cortex (PcunPCC).
In the next step, the spatial grouped-average EFs as an indicator of modulation intensity was calculated in all of the networks and inter-individual variabilities across the subjects were quantified by standard deviations (SDs). To exclude the networks with low average tDCS-induced EFs from the further analyses in the next steps, we defined a threshold. 99th percentile of the EF (an indicator of hotspot inside the brain) was calculated for each subject and, inspired by 19 , 50% of the lowest 99th percentile among participants was defined as the threshold. Subsequent calculations and comparisons were only performed at networks with mean tDCS EF magnitude higher than this threshold.
For statistical analysis, t-tests were used to examine significant differences between the two montages. Averaged EFs in each parcellated regions were compared between the two arrangements. False discovery rate (FDR) correction was used to correct P-values as necessary to overcome the multiple comparisons problem and P-value < 0.05 was considered to be significant. All data reported as mean ± standard deviation (SD).

Results
Surface-based CHMs were simulated for 66 participants with F4-Fp1 and F4-F3 montages. Personalized CHMs for F4-Fp1 and F4-F3 montages over all subjects are available in Fig. 2 which indicates a visible variation in tDCS induced EFs within a montage among participants. The spatial global maximum was calculated in the whole brain for each individual as an indicator of the EF hotspot. As shown in Fig. 3, a considerable variability for maximum EF is found among subjects but there is no significant difference between F4-Fp1 (0.41 ± 0.07 V/m; range from 0.28 to 0.58) and F4-F3 (0.39 ± 0.07 V/m; range from 0.25 to 0.59) montages in terms of global maximum EF.
For group-level analysis, the individual EFs were transformed into the standard fsaverage template. Transformation to standard space underestimated the global peak on average by 0.3%, showing that peaks are smoothed  www.nature.com/scientificreports/ a little in the transformation to standard space. Group-averaged EF among 66 participants in standard space can be found in Fig. 1 part (E) for both montages and it can be found that the highest EF intensity occurs in the frontal lobe.
Results at the network level. Figure 4 represents the amount of mean EF intensity at seven large-scale brain networks in the right and left hemispheres. The exact amount of mean EF intensity and P values for each network are available in Table S1 in Supplementary Materials. The results show that, regardless of montage choice here, the highest EF is produced in limbic system (F4-Fp1: left = 0.1955 ± 0.04, right = 0.155 ± 0.03; F4-F3: left = 0.1540 ± 0.03, right: 0.1544 ± 0.03) and lowest EF can be found in visual network (F4-Fp1: left = 0.0523 ± 0.01, right = 0.0532 ± 0.01; F4-F3: left = 0.0473 ± 0.01, right = 0.0488 ± 0.01) in both montages. In the left hemisphere, differences between two montages were statistically significant in Vis (P Corrected < 0.01), Dor-sAttn (P Corrected < 0.01), and Limbic (P Corrected < 0.001) networks. In the right hemisphere, the difference between the two montages was significant only in Vis (P Corrected < 0.01) network. Based on the predefined threshold, as shown with a green horizontal line in Fig. 4, tDCS induced EFs in Vis, SomMot, and DorsAttn networks are not above the threshold. Therefore, the left Limbic is the only network that has a significant difference between the two montages and EF intensity crosses the threshold in this network. Inspired by 19 , we considered 50% of the lowest peak as a threshold (Th = 0.125 V/m) to exclude the networks with low average EFs from further analyses in the next steps. The threshold is now represented by the horizontal green line in all of the bar graphs. Inspired by previous group-level CHM studies 19,31,32,[48][49][50][51] , we transformed all of the personalized CHMs to the standard space. Distortion of the EFs by transformation to the standard space can be a challenge in the group-level analysis of CHMs. To make sure that transformation did not highly affect our results we replicated our pipeline in the subject-space (using atlas-based parcellation in the subject-space instead of extracting EFs F4-F3 (red bars) electrode montages in the left (first row) and right hemisphere (second row). Labels below the horizontal axis determine the name of each network in Yeo7 parcellation and small brains next to the labels represent each network in fsaverage space. Significant differences between the two montages are shown above the bars based on the t-test with FDR correction threshold at P < 0.05. The horizontal green line indicates EF threshold (50% minimum value of 99th percentile of EFs in whole-brain analysis across the subjects). EF electric field, SD standard deviation, FDR false discovery rate, Vis visual network, SomMot somatomotor network, DorsAtn dorsal attention network, VentAtn ventral attention network, Limbic limbic system network, ECN executive control network, DMN default mode network. www.nature.com/scientificreports/ from standard space) for a part of our results. We compared our subject-space results with standard-space and we found no significant differences between these two methods at the group-level (Fig. S1). We also checked the results subject by subject and we found that there are very slight differences between EFs extracted from ROIs in the subject space compared to the standard space ( Fig. S2) without reaching to any statistical significance threshold.
The extent of inter-individual variability in the EF was calculated by average induced EF in each network that can be assumed to be normally distributed (Shapiro-Wilk test P > 0.05). Figure 5 shows the corresponding histograms for F4-Fp1 (blue) and F4-F3 (red) montages over 66 participants in the left and right hemispheres. Histograms display inter-individual variability in each subject's mean EF in the targeted network. We present the distribution of the EF strength in the networks above the threshold and results showed that inter-individual variance is relatively high in all of the networks. Except for the left Limbic, for all other networks, we found a relatively similar distribution in both montages.
Furthermore, inter-individual differences between subjects in the range of EF alterations by changing the electrode montage are represented in Fig. 6.
Results at the sub-network level. According to the results of the network analysis in the previous step, we showed that applied EFs in the four networks are above the threshold value. In Fig. 7, average EFs intensities in the main nodes of these networks are represented separately for each hemisphere. The numerical values of averaged EFs in each node can be found in Table S2 in the Supplementary Materials. As indicated in Fig. 7, in most areas, the EF intensity generated by F4-Fp1 montage is higher than F4-F3. Statistical results indicate that there is no significant difference between the two montages in the nodes with EF above the threshold in the right hemisphere. Nonetheless, in the left hemisphere there are some main nodes with significant differences between two montages including FrOperIns (F4-Fp1: 0.1514 ± 0.03, F4-F3: 0.1385 ± 0.03; P < 0.01) node in Ven-tAttn, OFC (F4-Fp1: 0.2474 ± 0.05, F4-F3: 0.1953 ± 0.04; P < 0.001) and TemPole (F4-Fp1: 0.1342 ± 0.03, F4-F3: 0.1073 ± 0.02; P < 0.001) in Limbic, and ECN-PFC (F4-Fp1: 0.1937 ± 0.04, F4-F3: 0.1661 ± 0.04; P < 0.001) that received EFs above the threshold values.

Discussion
In this study, we modeled the network selectivity of two tDCS protocols (F4-Fp1 "unilateral" and F4-F3 "bilateral" montages) that are commonly used for stimulating DLPFC. We tested the individual head models for the electric field (EF) at the group level to examine which parts of the large-scale brain networks are most likely to be stimulated by considering two conventional electrode montages. We showed that the placement of the cathode electrode can significantly change EF distribution patterns in specific nodes in large-scale networks even if there is no significant difference in the whole brain network. We found that in the F4-Fp1 montage, the limbic system network in the left hemisphere receives a significantly higher electrical dose. Within the network nodes, frontal-operculum-insula nodes in the ventral attention network, orbitofrontal cortex and temporal pole in the limbic system, and DLPFC node in the executive control network significantly receive higher EFs by F4-Fp1 arrangement. Figure 2 visualizes the inter-individual variability in the EFs. A wide range of variation in the EF distribution patterns over 66 participants supports the previous findings on the importance of considering the individual computational head modeling (CHM) 18,32,38,51 . Individualized head models for a group of subjects with www.nature.com/scientificreports/ surface-based registration to standard template allowed us to use standard functional atlas parcellation for network-level analysis of the head models.
In this work, we extracted the EF intensity inside the large-scale brain networks and their nodes (Figs. 4, 7). Calculation of EF intensity in brain regions requires a decision about which EF measures (e.g. maximum, average, median) are appropriate to use. Gomez et al. by using network parcellation of the cerebellum reported that in contrast to maximum, functional network analysis becomes important when average EF is used 38 .
Our results conform with previous computational modeling studies that cathode location has a significant impact on EF distribution patterns 49,52 . Previous findings reported that the effects of tDCS depend on the locations of both active and reference electrodes 28,30 .
Our results indicate that among the networks with EF intensity above the defined threshold, the differences between the two montages were significant only in the left limbic system. Cathode locations in both arrangements are spatially close but nonetheless cathode over the left orbit is closer to the orbitofrontal cortex (OFC) than cathode over left DLPFC (F3). Therefore, F4-Fp1 induced higher EFs than F4-F3 inside this region as the main node of the limbic system. Besides between montage differences, as shown in Fig. 4, we also observed within montage variance in terms of averaged EFs inside a network (Figs. 5, 6). For instance, by changing the electrode location from the left supraorbital area to F3 according to the 10-20 EEG electrodes placement standard, there were some subjects whose EFs changed extremely little or oppositely compared to the grouped-average results. For example, changes of EFs in the left limbic system, that showed a significant difference between two montages, range from − 2 × 10 -5 to 6.9 × 10 -2 [V/m] and there is one person (1.5% of total participants) whose EFs changed in opposite direction compared to the grouped-average EFs. However, EF changes in the right limbic system range from − 2.8 × 10 -2 to 1.7 × 10 -2 , and there are 19 subjects (28.8% of total participants) whose EFs changed in opposite direction compared to the grouped-average. With respect to previous simulation studies, these inter-individual variations in EF distribution patterns could be explained by different anatomical factors that can change the actual current received by different parts of the cerebral cortex 14,44,51 . The current distribution of the tDCS is influenced by microarchitectural features of the brain and variation in the cortical anatomy such as gyrus and sulcus patterns can change induced EFs between individuals 53,54 .
EF calculations in the main nodes are depicted in Fig. 7 and numerical results can be found in Table S2 in Supplementary Materials. According to the brain's compensatory response to intervention through the neural networks, a simple hypothesis is that modulating multiple nodes of a network would be more effective than modulating a single node in that network 55 . As shown in Fig. 7, the intensity of the induced EFs in distinct nodes of a network is different. Grouped-average EF intensity is weak in some nodes in a network since they www.nature.com/scientificreports/ are located farther from stimulating electrodes. For example, DLPFC, as the main node of the executive control network (ECN) was targeted in both montages in our simulations. This node was named ECN-PFC in Fig. 7 and received the highest EF inside the network. Statistical results indicated that induced EF by F4-Fp1 in DLPFC was significantly higher than F4-F3. However, EF intensity in other parts of the ECN, precuneus-cingulate, and temporal-parietal nodes, is considerably low. This opens potentials for multifocal tDCS montages stimulating all parts of a specific distributed network 37 . We also found that the PFC node in the DMN network received the highest EFs inside the networks and the amount of EFs in other nodes of this network is considerably low. The results from a tDCS-fMRI study conducted Labels below the horizontal axis denote the name of main nodes inside the largescale network and small brains next to the labels represent the spatial location of the main nodes in fsaverage space. Significant differences between the two montages are shown above the bars based on the t-test with FDR correction threshold at P < 0.05. The horizontal green line indicates EF threshold (50% minimum value of the 99th percentile of EFs in whole-brain analysis across the subjects). EF electric field, SD standard deviation, FDR false discovery rate, VentAtn ventral attention network,  Fig. 7, our calculations showed that temporoparietal and PCC nodes in DMN received a small amount of EFs. Perhaps a significant diminishing in functional connectivity of these nodes reported by Shahbbie et al. is related to low EF intensity in these nodes that have not yet been investigated. Other than ECN-PFC, analysis of the main nodes revealed that, although the left limbic system is the only network that shows a significant difference between two montages, however, frontal-operculum-insula node in the ventral attention network and orbitofrontal cortex and temporal pole nodes in the limbic system all in the left hemisphere are significantly different between these two montages. These significant differences may become a matter based on the stimulation goal. For example, the frontal operculum in the ventral attention network is a key node in the cognitive control system. In a combined TMS-fMRI study, it was reported that the frontal operculum regulated the level of activity in the occipitotemporal cortex. By applying TMS, they found that stimulation of this region decreased top-down selective attentional modulation in the occipitotemporal cortex 57 . Based on our calculations, F4-Fp1 generated significantly higher EF inside this node of the ventral attention network. Therefore, it might be a more effective montage than F4-F3 for stimulating frontal-operculum-insula which is important in human cognitive tasks including directed attention tasks. However, this is an empirical question for future trials.
As shown in Fig. 7, both OFC (within the limbic network) and DLPFC (PFC node in ECN) were stimulated strongly regardless of electrode montage and greater EFs induced by the F4-Fp1 arrangement. These two nodes play important role in different cognitive functions such as decision-making behavior and inhibitory control. A previous tDCS study compared the effects of unilateral (i.e., F4-Fp1 or F3-Fp2) and bilateral (i.e., F4-F3 or F3-F4) DLPFC stimulation on risk-taking and decision-making behaviors. Similar to our simulations, they used 5 × 7 large electrode pads and the current strength was 2 mA. They reported that based on Balloon Analog Risk Task (BART) as a measure of risk-taking propensity, bilateral stimulation of the DLPFC (anode-cathode over F4-F3 or F3-F4), compared to sham, significantly reduced risk-taking behavior. While no statistically significant difference in risk-taking score was discovered between sham and unilateral DLPFC 43 . In a tDCS-fMRI study neural effects of bilateral tDCS (F3-F4) on risk-taking behavior during BART task were examined. The authors reported that whole-brain connectivity of right DLPFC (under anodal electrode) correlated negatively with risk-taking aversion. It was also reported that OFC activity is negatively correlated with risk and pop rates in the BART task 58 . However, the authors have not investigated the relationships between the induced EFs by unilateral or bilateral montages and neural/behavioral outcomes. Testing the effect of higher levels of EFs in ECN and the limbic system generated by unilateral stimulation compared to the bilateral montage, that we found in our simulations, in relationship with neural and behavioral outcomes might give us a better understanding of the underlying mechanism of actions.
Previous research also showed that bilateral stimulation of DLPFC can affect the brain in a different way compared to unilateral stimulation during the temporal discounting task. The temporal discounting task is a measure of risky decision making, subjective value judgments, and the ability to delay gratification. No change in temporal delay discounting with bilateral tDCS over DLPFC was described by 59 . However, the unilateral DLPFC stimulation (F3-Fp2) modulated temporal discounting rate 60 . Future tDCS studies on cognitive functions critically depend on OFC and DLPFC activities like emotion, motivation, reward or valence processing 61 may like to consider the effect of inter-montage and intra-individual variations in induced EFs inside these nodes and their corresponding network on neural and behavioral outcomes based on the analysis pipeline and the initial findings we have introduced in this study.

Limitations and future direction
Our study has some limitations that could be addressed in future research. The main limitation is that our results are purely based on simulation. We did not consider other outcome measures such as neural activity/connectivity or behavioral outcomes. We considered averaged EF as an indicator of stimulation intensity. Open questions remain about how to relate EFs with behavioral or neurophysiological changes. The relationship between EF intensity and neural/behavioral responses is not a trivial matter and there is still no clear understanding of the underlying tDCS mechanism of action at the network level 62 . Integrating EF modeling and neuroimaging information such as functional connectivity that incorporates the interactions within and between large-scale functional networks can provide for a better understanding of tDCS effects in terms of dose-response relationships. For example, Esmaeilpour et al. reported a significant correlation between EF obtained from CHMs and changes in fMRI signals 62 . Kasten et al. revealed a direct link between EF variability and alpha band frequency changes as a proxy indicator of neural activity in a transcranial alternative current stimulation (tACS) study combined with MEG 23 . Combining CHMs with neuroimaging data provide critical insights about tDCS mechanistic effects. However, it is still in an early stage of development and more analysis is needed to increase our understanding of induced EFs and tDCS outcomes in general.
Furthermore, we did not consider network interaction in our analysis pipeline. To our best knowledge, there is no published attempt to find the relationships between induced EFs and interaction between networks that have a critical role in cognitive processes yet. For instance, previous tDCS studies in the treatment of addictive disorders showed potential effects by targeting networks known to be dysregulated in SUDs such as ECN for executive control and DMN for internal ruminations 63 . These two functional networks typically act in counterbalanced order. The hypothesis is that by applying anodal stimulation to DLPFC the level of activity/connectivity should be increased in ECN while it should be decreased in DMN 64  www.nature.com/scientificreports/ (F4-F3) decreases functional connectivity in DMN especially in posterior parts of the network including the middle temporal gyrus, right precuneus, and right PCC 56 . As represented in Fig. 7, our calculations showed that temporoparietal and PCC nodes in DMN receive a small amount of EFs in F4-F3 montage. Integrating CHMs with functional data at the network level can shed light on how large-scale functional networks interact with each other by applying tDCS. Although we calculated EFs in the main nodes of the networks, it is remained unclear which nodes outside the DLPFC are actually important for enhancing tDCS effects on intervention outcomes. The contribution of EF intensity in each part of the network with tDCS outcomes is ambiguous and there is not enough evidence to conclude which nodes should be stimulated, and which ones should be inhibited to enhance stimulation outcome measures. More physiological and behavioral information is needed to find a relationship between the main components of the networks and tDCS inter-individual responses. This unresolved question would be critical, especially for optimal dose selection and dose customization.
Furthermore, in this study, we only simulated conventional large electrode pads to determine induced EFs at the network-level. Focal stimulation of the network nodes in a group of subjects by using conventional electrodes is difficult. Our study suggests that the probability of interaction between large-scale networks should be considered more carefully using conventional tDCS. Stimulating many nodes of different networks because of the diffusivity of the current makes complexities in modeling the network level contribution in the conventional tDCS montages. Using protocols that produce more focal stimulation such as high-definition (HD) 24 , or multiarray electrodes could potentially provide the possibility to control modulatory effects at the network level 37 . In this context, simulation of the HD and multi-array electrodes might bring new insights for the network-level analysis of CHMs.
Another limitation is that our simulations are focused on EF intensity (magnitude of the EF) inside the networks, which is informative of the EF strength. In the current study, we ignored EF components radial to the cortical surface (normal components of the EF), which is informative of EF direction. The normal component of the EF, which is perpendicular to the cortical surface, reflects injected current either entering or leaving the cortex. Previous researches suggest that this component causes polarity-specific effects including facilitatory (anodal like) or inhibitory (cathodal like) effects 5 .
Atlas-based parcellation provides information from a group of subjects that may be different from the person's resting-state functional networks. Parcellation of CHMs based on the individual's functional connectivity might be more accurate to determine EFs inside the main nodes of large-scale-networks. Furthermore, we have only focused on the organization of large-scale distributed networks obtained from the atlas of resting-state (task-independent) functional networks. Calculation of the EFs according to task-evoked connectivity was not considered in the present study. Investigation of the induced EFs over task-based networks may help to target central brain regions that play pivotal roles in the task performance. For example, a previous study showed that under anodal tDCS (F7-Fp2 montage), during a word generation task, changing connectivity in a given taskrelated network provides the basis for enhanced neural efficiency in highly specific brain areas critical for task performance. These areas might be located under the stimulation site or distant connected regions 65 . These findings suggest that induced EFs in the major hubs of the task-related network may help understand the mechanistic effects of tDCS. Moreover, targeting main hubs in the task-related networks and functionally connected regions using CHMs may increase the effectiveness of tDCS during task performance since it is hypothesized that tDCS acts on the networks and brain regions concurrently active.
EF modeling in this study is performed based on many assumptions. We used previously established isotropic conductivities for tissues, as is common in computational studies. Suh et al. reported that the anisotropic skull and WM conductivities significantly affect EF distribution patterns 66 . However, Shahid et al. revealed that anisotropy does not modulate significant changes in comparisons across montages 67 . As mentioned in another study, modeling anisotropy is important when considering deeper target regions inside the brain 68 . Furthermore, our tissue conductivities were completely static, and we did not consider the effects of applying current on tissue conductivities. As mentioned in 42 , electrodes, tissue, and their interfaces are not merely resistive and can be changed by the applied current intensity that may vary over time. However, these nonlinear effects have not yet been considered in tDCS modeling studies so far.
Next, we used automatic tissue segmentation and outcomes are strongly dependent on the quality and resolution of T1 and T2-weighted images. However, the quality of segmentation results was evaluated visually sliceby-slice to ensure correct tissue classifications. Furthermore, warping from native space to standard fs average space may cause a certain degree of information loss, but it's necessary for group-level analysis.
Future work is needed to determine whether induced EFs in the main nodes of the large-scale brain networks link to changes in neural activities. It should be investigated that when a single brain region is targeted for stimulation, how other nodes of that region's network and main components of the other networks interact with the targeted area. Combining CHMs with neuroimaging and behavioral data provides potential insights for understanding how induced EFs interact with the functional activity of the brain at the network level and when stimulation becomes more effective based on current density inside the main components of the networks. Our simulation results from a group of participants with MUDs could be integrated with their clinical and physiological data to determine interindividual variations. Additionally, it would be possible to propose a customized montage by considering the initial state of each individual (e.g. functional connectivity or behavioral measurements at the baseline) and distribution of the current in personalized CHMs. Such effects have not yet been incorporated into previous brain stimulation modeling studies.

Conclusion
In this study, we showed that how current flows through the large-scale networks based on group-level analysis of CHMs. In addition to inter-individual variability, we showed that the spread of the EFs in main nodes of the large-scale networks was significantly different between two similar electrode montages for DLPFC stimulation and cathode location can change EF distribution patterns at the network-level especially in the limbic system. The calculation of EFs in the networks and main nodes would be informative since brain stimulation is a network phenomenon. Knowledge about functional connectivity within and between large-scale networks and the EF distribution inside the brain networks might help to resolve ambiguity about how much tDCS effects spread across the brain. The proposed method in this study suggests that a network parcellation of CHMs at the grouplevel can be used in future studies to understand how tDCS affects large-scale networks and how the results might vary depending on inter-individual variations and electrode arrangements.

Data availability
The raw database and CHMs generated for this study are available on request to the corresponding author.

Code availability
EF intensity results and MATLAB codes can be found here https ://osf.io/3ukbc /.