Addressing the sources of inter-subject variability in E-field parameters in anodal tDCS stimulation over motor cortical network

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ This is the Accepted Manuscript version of an article accepted for publication in Physics in Medicine & Biology. IOP Publishing Ltd is not responsible for any errors or omissions in this version of the manuscript or any version derived from it. The Version of Record is available online at https://doi.org/10.1088/1361-6560/ad5bb9 Accepted Manuscript under a CC BY-NC-ND licence. 1 Pablo Franco-Rosado1,2,3,***, M. Amparo Callejón3,4, Javier Reina-Tosina3, Laura M. Roa3, Juan F Martin-Rodriguez1,2,5, Pablo Mir1,2,5 1. Unidad de Trastornos del Movimiento, Servicio de Neurología, Instituto de Biomedicina de Sevilla, Hospital Universitario Virgen del Rocío/CSIC/Universidad de Sevilla, Sevilla, Spain. 2. Centro de Investigación Biomédica en Red sobre Enfermedades Neurodegenerativas (CIBERNED), ISCIII, Madrid, Spain. 3. Grupo de Ingeniería Biomédica, Escuela Técnica Superior de Ingeniería, Sevilla, Spain. 4. Servicio de Otorrinolaringología, Hospital Universitario Virgen Macarena, Sevilla, Spain 5. Departamento de Psicología Experimental, Universidad de Sevilla, Sevilla, Spain 6. Departamento de Medicina, Facultad de Medicina, Universidad de Sevilla, Sevilla, Spain ***: Corresponding author (pfranco[email protected]) 2 Addressing the sources of inter-subject variability in E-field parameters in anodal tDCS stimulation over motor cortical network ABSTRACT Background: Although tDCS constitutes a non-invasive neuromodulation technique with promising results in a great variety of applications, its clinical implementation is compromised by the high inter-subject variability reported. This study aims to analyze the inter-subject variability in electric fields (E-fields) over regions of the cortical motor network under two electrode montages: the classical C3Fp2 and an alternative P3F3, which confines more the E- field over this region. Methods: Computational models of the head of 98 healthy subjects were developed to simulate the E-field under both montages. E-field parameters such as magnitude, focality and orientation were calculated over three regions of interest (ROI): M1S1, SMA and preSMA. The role of anatomical characteristics as a source of inter-subject variability on E-field parameters and individualized stimulation intensity were addressed using linear mixed-effect models. Results: P3F3 showed a more confined E-field distribution over M1S1 than C3Fp2; the latter elicited higher E-fields over supplementary motor areas. Both montages showed high intersubject variability, especially for the normal component over C3Fp2. Skin, bone and CSF ROI volumes showed a negative association with E-field magnitude irrespective of montage. Grey matter volume and montage were the main sources of variability for focality. The curvature of gyri was found to be significantly associated with the variability of normal E-fields. Conclusions: Computational modelling proves useful in the assessment of E-field variability. Our simulations predict significant differences in E-field magnitude and focality for C3Fp2 and P3F3. However, anatomical characteristics were also found to be significant sources of E-field variability irrespective of electrode montage. The normal E-field component better captured the individual variability and low rate of responder subjects observed in experimental studies 1. INTRODUCTION Transcranial direct current stimulation (tDCS) is a non-invasive neurostimulation technique aimed at modulating cortical excitability through the application of low-intensity (1-2 mA) electric currents into the brain through two or more electrodes placed on the scalp [1], [2], [3]. tDCS has shown promising results in its application to a large variety of neurological disorders involving cognitive and motor functions [4], [5], [6], [7]. However, open questions remain, such as the delineation of the true targeted brain regions responsible for the tDCS-induced 3 physiological effects [8], [9], [10] and the high inter-subject variability observed in experimental studies [11], [12], [13]. Specifically, when stimulating the motor cortex in tDCS, different electrode montages have been proposed. For instance, in addition to the classical montage proposed by Nitsche et al. [1], [14] with the anode over the left motor cortex and the cathode in the contralateral supraorbital region; an alternative parieto-frontal montage for which the E-field distribution is expected to be mostly confined ipsilateral over the targeted motor area has also been tested [15], [16]. Despite the differences observed in E-field distribution for these two montages, both have shown to successfully modulate cortical excitability in the motor cortex at a group level [1], [17], [18]. However, experimental studies have also reported significant inter-individual variability and a high rate of non-responder subjects [19], [20]. In addition, regarding tDCS targeted regions, experimental studies have suggested that, in addition to the motor cortex, the supplementary motor cortex (SMC) may play a crucial role in various pathologies such as in Parkinson disease, where a decrease in activity is observed in comparison to healthy subjects, resulting in altered gait, sequential performance and temporal processing [21]. Preliminary studies indicate that the stimulation of SMC can improve Parkinson's symptoms [22], [23], [24], thus highlighting the importance of evaluating the impact of neuromodulation in these supplementary regions for the effective design of therapeutic options. Realistic computational modelling approaches have proved essential tools to increase our understanding of the biophysical basis of tDCS by predicting the magnitude and distribution of the induced electric field (E-field) in the brain [25], [26], [27]. For instance, the use of realistic models has provided a path for the study of the influence of individual anatomy on E-field distribution as a major cause of the experimental variability observed [18], [28], [29], [30]. They have also put the conventional criteria for electrode placement into doubt when predicting considerable E-field spread over brain areas relatively far from the presumed targeted regions [31], [32], [33], [34]. However, the relation between the predicted E-fields and the induced neurophysiological effects remains under discussion, since they do not completely explain the divergence found in different experimental studies [35], [36], [37]. Specifically, it is not completely understood how and to what extent individual differences influence E-field parameters irrespective of electrode montage, or conversely, whether specific electrode montages can be used to minimize inter-subject variability [38], [39]. 4 In this work, we develop a full-head modelling framework with the objective to analyze and compare relevant E-field parameters, such as magnitude, focality, and orientation (through tangential and normal E-field components), under the two electrode montages mentioned above: the classical C3Fp2 and an alternative ipsilateral P3F3 montage. In addition, we expand the analysis of the E-field [40], [41], [42] over the sensorimotor region M1-S1 to cover the supplementary motor cortex (SMC) composed of the secondary motor areas SMA and preSMA. We also predicted individual supra-threshold stimulation intensities based on these full-head E- field models. In resume, the main issues and research questions investigated in this study are: 1) analysis of inter-individual differences in predicted E-fields for each electrode montage: does the alternative montage P3F3 reduce the variability of the simulated E-field?; 2) does the alternative montage lead to optimized E-field parameters in terms of magnitude and focality over the main targeted M1-S1 and over the secondary SMA and preSMA regions?, and 3) which are the main anatomical characteristics responsible for the variability observed in E-fields and individual stimulation intensities?. 2. MATERIAL AND METHODS This section describes the material and methods followed to develop the modelling pipeline framework shown in Fig. 1, structured as follows: Section 2.1 describes the characteristics of the full head models obtained and Section 2.2 the two electrode montages emulated. Section 2.3 exposes the details of the simulations conducted using the finite-element method. Section 2.4 details the characteristics and parameters of the E-fields simulated and Section 2.5 the statistical analysis performed. 2.1 Anatomical Head Models Head models were created for 98 healthy volunteer subjects (52 women, mean age: 32 ± 8 years) from T1 (TE = 3.8 ms, flip angle = 8◦, 1 x 1 x 1 mm voxel) and T2-weigthed (TE = 252.8 ms, TR=2500 ms, flip angle=90◦, 1 x 1 x 1 mm voxel) MRIs acquired with a Philips Ingenia 3T scanner. All recruited volunteers were informed of the purpose of the study and signed the written informed consent form. The use of the same MRI sequence in all subjects avoids a source of variability in the construction of the full head models based on tissue segmentation. Models were constructed using the automated headreco workflow available on SimNIBS 3.2 software package (https://simnibs.github.io/simnibs/build/html/index.html) [43]. The use of high-resolution meshes allows an accurate characterization of tissue topology and a detailed analysis of the induced E-field over cortical gyri and sulci [44], [45]. In our models, final head 5 meshes contained more than 19 million tetrahedrons on average, with a mean edge size of 1.22 mm, assigned to different head tissues including scalp, bone, vitreous humor (eyes), air cavities, cerebrospinal fluid (CSF), grey matter (GM) and white matter (WM). 2.2 Electrode Montages Electrodes were modelled as rectangular rubber electrodes with an area of 4.5 x 6.5 cm² (DC- Stimulator, Eldith GmbH) and a thickness of 1 mm, respectively built on a 3mm-thick layer of conductive gel. A metal connector was modelled 1 mm below the electrode with an area of 1.6 x 2.3 cm². Two electrode configurations for motor cortex stimulation were considered: conventional C3Fp2 montage, with anode and cathode placed over C3 and Fp2, respectively [13], [46]; and P3F3 [16], an alternative posterior-anterior configuration with both electrodes approximately equidistant to the precentral gyrus, with the anode on P3 and the cathode on F3. Electrode positions refer to the 10/10 EEG system. We selected these two electrode montages due to their expected differences in E-field magnitude, distribution and orientation over the targeted motor region [57], presumably allowing us to comprehensively analyze and compare these two main relevant E-field parameters: larger E-field magnitude is expected between electrodes under P3F3 montage rather than directly below them as in C3Fp2 [47], [48]. However, the orientation of the E-field is expected to be more tangential in P3F3 compared with C3Fp2 [49]. Therefore, these two electrode montages serve us as representative examples to guide the research questions addressed in this study. Figure 1. Modelling pipeline for tDCS simulation and analysis. Steps 1-3 were carried out in an automatized way using SimNIBS 3.2 software: (1) segmentation of T1 and T2-w MRIs and generation of volume conductor head models, (2) placement of electrodes under two different configurations, C3Fp2 and P3F3, and (3) calculation of the E-field using FEM. (4) Group and individual analysis of results: (4a) the simulated E-fields were exported to the normalized fsaverage surface to calculate spatial correlation across subjects 6 for the norm of the E-field and for their tangential and normal components. (4b) Individual analysis in the GM volume in three regions of interest, M1-S1, SMA, and preSMA. Finally, (5) statistical analysis and comparison of magnitude and focality of the E-field for each ROI and electrode montage is performed. 2.3 FEM Calculation SimNIBS was used to solve the physics for tDCS under quasi-stationary conditions, i.e., Laplace equation (𝛻𝛻 ∙ (𝜎𝜎𝛻𝛻𝜎𝜎)= 0, to obtain the E-field (𝐸𝐸 �  =−𝛻𝛻𝜎𝜎) and the current density (𝐽𝐽 = 𝜎𝜎𝐸𝐸 �  ) in each element domain. A stimulation intensity of 1.5 mA was selected, as in previous studies [50], [51]. Tissues were considered as isotropic and homogeneous with conductivity values available by default in SimNIBS: 0.126 S/m for WM, 0.275 S/m for GM, 1.654 S/m for CSF, 0.01 S/m for bone, 0.465 S/m for scalp [52], 0.5 S/m for eyes [33], 1 S/m for the conductive gel layer [53], and 29.4 S/m for the rubber electrodes (NeuroConn electrodes: Wacker Elastosil R 570/60 RUSS). 2.4 E-field simulation 2.4.1 Normalized fsaverage surface In order to analyze and compare the spatial distribution of the E-field for C3Fp2 and P3F3 montages, the E-fields obtained for each head model were exported onto the normalized surface fsaverage, which is a high-resolution triangular mesh of 327k nodes representing the topology of the cortical surface [54], [55]. The average E-field magnitude (calculated over 98 subjects) was computed at each surface node. Over this surface, the multimodal cortical parcellation of the Human Connectome Project (HCP) was considered to delineate the sensorimotor region through the connection of regions 6 to 9, as described in [56]. 2.4.2 Individual ROIs A modelling workflow was specifically designed to extract and analyze different ROIs in the individual GM volume models, including the sensorimotor area M1-S1, the supplementary motor area (SMA), and the preSMA. First, a binary mask was obtained for the M1-S1 region using the image calculator fslmath, through the combination of the precentral and postcentral gyrus of the left hemisphere of the MNI152 Harvard-Oxford cortical structural atlas (https://cma.mgh.harvard.edu/) and the mri_binarize module available on FSL v.6.0 software (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) [57]. The binary masks corresponding to the SMA and preSMA regions were extracted from the Sallet’s dorsal frontal connectivity-based template [58], also available in FSL. The resulting masks were manually edited to ensure that the corresponding tetrahedral meshes did not overlap. The resulting masks were transformed from MNI to the subject’s own space and subsequently converted to triangular meshes. Finally, 7 through customized MATLAB-based scripts (R2019a), the intersection of these triangular meshes with the GM volume of each head model was obtained. 2.4.3 Individual supra-threshold stimulation intensities Following the procedure reported in [59], we assumed E-field values larger than 0.15 V/m as thresholds values for neuromodulatory effects over the hand knob, based on previous modelling works that reproduced clinical studies where these effects were observed [60], [61]. As such, we calculated the individual stimulation intensities which led to E-field values larger than 0.15 V/m in the hand-knob region for each of the two montages selected. For volume-level analysis, the hand-knob region was defined as the left A4ul (upper limb) region of the Brainnetome atlas [62], which correlates functional connectivity and anatomical information. The same procedure described in section 2.4.2 was applied to determine the volumetric mesh of this region. The description of the hand-knob region for the surface-level analysis was done by combining areas 6mp, 6d, and 4 from the HCP. Finally, for the sake of comparison with experimental studies, in which the most frequent range for anodal stimulation ranges between 0.5 and 2 mA [63], [64], usually constrained due to safety recommendations [65], [66]; we compared the proportion of subjects who were able to reach the 0.15 V/m E-field neuromodulatory threshold with anodal intensities lower than or equal to 2 mA for the two electrode montages simulated. This variable was defined as the responder rate (RR) for each electrode montage. 2.5 Data analysis 2.5.1 Analysis of spatial correlation across subjects at the fsaverage surface To asses inter-subject variability, the spatial correlation of normal and tangential E-fields over the sensorimotor region on the fsaverage surface was calculated. The mean correlation and standard deviation (SD) obtained for the 98 subjects were compared between the two electrode montages. The coefficient of variation, CV (%), was used to quantify the inter-subject variability or deviation with respect to the mean of the group. 2.5.2 Analysis of E-field parameters at individual ROIs Second, individual mean and maximum E-field values for each ROI and montage were calculated. In order to avoid spurious values, the maximum E-field was filtered by the 99th percentile. In the literature, focality is usually defined as the relative volume inside the ROI exposed to an E-field strength higher than a particular percentile value [16], [67], [68]. In this study we defined a percentile value of 70 to analyze the focality of the E-field in the individual 8 ROIs. The normality of individual E-field parameters was assessed through Lilliefors test. Statistical differences between the two electrode montages were analyzed by using paired Student’s t-test or Wilcoxon signed rank test depending on distribution. Effect sizes measures were reported through Cohen’s d value [69]. 2.5.3 Analysis of the effect of individual anatomical characteristics on E-field parameters and stimulation intensity To analyze and quantify the effect of anatomical characteristics on predicted E-fields irrespective of electrode montage, linear mixed models were calculated, considering electrode montages C3Fp2 and P3F3 as level-one units nested in subjects (level-two units). The dependent variables considered were average and maximum E-field and focality in the targeted M1-S1 region. Anatomical characteristics such as total intracranial volume (TIV), skin and skull volumes below the stimulating electrode; as well as CSF, GM, and WM volumes in the M1-S1 region were considered as covariates, as well as sex and age. Beta coefficients, confidence intervals (Cis) and corresponding p-values were reported. Coefficient of multiple determination value (R2) was used to assess model fitting. In addition, semi-partial R2 were calculated as measures of effect size for fixed effects. A p-value of 0.001 was considered significant in our study to detect more robust statistical associations. Additionally, in order to quantify the contribution of covariates to the average deviation in predicted E-field parameters separately for each electrode montage, simple and multiple linear stepwise regression models (backward and forward) were calculated following the deviance criterion with a threshold value of 0.001. With the same characteristics, the analysis was repeated with supra-threshold stimulation intensity as the dependent variable, assessing the effect of anatomical characteristics on its predicted magnitude. Finally, an analogous linear regression analysis was performed to evaluate the effect of morphological parameters of the sensorimotor HCP cluster for both the tangential and normal components of the predicted E-field and on the hand knob region for the individual suprathreshold stimulation intensities. The following morphological parameters were considered as covariates: cortical thickness in mm [70], gyrification (a measure which quantifies cortical folding complexity [71], fractal dimension (an index ranging between 0 and 4 to quantify cortical surface complexity), and sulcus depth (defined as the euclidean distance between the cortical surface and its convex hull) [72]. These parameters were calculated using the CAT12 (https://neurojena.github.io/cat/) segmentation included in headreco pipeline, which allowed us to obtain the distribution of these morphological values for each node of the fsaverage surface. The evaluation of such parameters describing cortical topology is intended to delve into the 15 To further analyze how much variance was explained by each independent variable, we calculated a model with age, sex, skin, CSF, GM, and WM volumes as the only covariates (Supplementary Table 3, check Supplementary Fig. 6-9 for linear fitting and scatter plots for the four dependent variables and each of the two electrode montages). The results are consistent with those previously obtained with the linear mixed models. For average and maximum E- fields, a significant negative association was found with skin, bone, and CSF ROI volumes for the two electrode montages C3Fp2 and P3F3, with R2 values of 0.40 and 0.51 found for the skin. In the case of focality, a significant positive association was found with CSF, GM, and WM volumes only for P3F3 montage, with the largest association being found for GM volume (R2 = 0.35). Multiple linear regression models confirmed the durable effect of such anatomical variables in the E-field magnitude across the two electrode montages analyzed, accounting for 63-88 % of the inter-subject variance (Supplementary Table. 4). For the focality, no association with anatomical variables was found for C3Fp2. For P3F3, individual CSF and GM ROI volumes explained 56% of the inter-individual deviation under this montage (Supplementary Table 5). 3.4.2 Effect on tangential and normal E-fields Linear mixed models were developed to study the causes of variability for tangential and normal E-fields, using individual morphological cortical surface parameters (sulcus depth, gyrification, cortical thickness and fractal dimension) as covariates, in addition to age, sex and TIV. Results obtained for both components are listed in Table 5. Significant positive strong associations were found in both components for both montage (p<0.001), while gyrification shows a clear tendency towards significance (p < 0.01) for C3Fp2, reaching significance for P3F3 (p < 0.001). No other significant associations were identified for other covariates. The results and representation of the study of the variance explained by independent variables are shown in Supplementary Table 6 and 7 and Supplementary Figures 10-11. Table 4. Results of linear mixed model analyses for dependent variables such as average and maximum E-field (V/m) and focality (dl). Significant beta coefficients and p-values are shown in bold. Fixed effects Average E-field strength (V/m) beta 95%-CI p R² Intercept 0.4313 (0.385, 0.477) <0.001 Sex (males) 0.0015 (-0.005 ,0.0078) 0.63739 0 Age (years) 1.39e-05 (-0.0003, 0.0004) 0.93207 0 Montage (ref: C3Fp2) F3P3 0.0349 (0.032, 0.038) <0.001 0.22 TIV (dl) 0.0005 (5.4558e-5, 0.0011) 0.030761 0.05 Skin (dl) -0.0441 (-0.049, -0.0387) <0.001 0.24 Skull (dl) -0.0542 (-0.064, -0.0444) <0.001 0.11 CSF (dl) GM (dl) -0.0358 0.0001 (-0.043, -0.0284) (-0.0112, 0.0115) <0.001 0.98241 0.09 0 16 WM (dl) -0.0082 (-0.0187, 0.0022) 0.11898 0 Random effects sigma subject 0.0089 0.0069-0.0114 Total R2 0.95 Fixed effects Maximum E-field strength (V/m) beta 95%-CI p R² Intercept 0.7988 (0.6930, 0.90442) <0.001 Sex (males) 0.0008 (-0.01380, 015479) 0.96398 0 Age (years) -1.7e-05 (-0.0008, 0.0007) 0.91094 0 Montage (ref: C3Fp2) F3P3 0.0537 (0.0459, 0.0615) <0.001 0.16 TIV (dl) 0.0012 (4.28e-05, 0.0024) 0.058513 0.01 Skin (dl) -0.0796 (-0.0921, -0.0672) <0.001 0.25 Skull (dl) -0.1043 (-0.1269, -0.0816) <0.001 0.13 CSF (dl) GM (dl) WM (dl) -0.0550 -0.0101 -0.0137 (-0.0721, -0.0378) (-0.0364, 0.0162) (-0.0378, 0.0104) <0.001 0.45073 0.26215 0.06 0 0 Random effects sigma subject 0.017426 0.011972-0.025365 Total R2 0.86 Fixed effects Focality (dl) beta 95%-CI p R² Intercept -0.1267 (-0.5876, 0.3343) 0.5883 Sex (males) 0.0239 (-0.0399, 0.0878) 0.4609 0.00 Age (years) -0.0011 (-0.0044, 0.0021) 0.4896 0.00 Montage (ref: C3Fp2) F3P3 0.0454 (-3.02e-05, 0.0909) 0.0502 0.02 TIV (dl) 0.0013 (-0.0040, 0.0066) 0.6239 0.00 Skin (dl) 0.0344 (-0.0201, 0.0888) 0.2144 0.01 Skull (dl) 0.0154 (-0.0835, 0.1143) 0.7598 0.00 CSF (dl) GM (dl) WM (dl) -0.0732 0.2158 -0.0255 (-0.1480, 0.0016) (0.1012, 0.3304) (-0.1310, 0.0796) 0.0545 <0.001 0.6327 0.01 0.05 0.00 Random effects sigma subject 2.41e-07 - Total R2 0.24 Table 5. Results of linear mixed model analyses for morphological parametres of cortical surface: depth (mm), thickness (mm), gyrification and fractal dimesion. Significant beta coefficients and p-values are shown in bold. Fixed effects Tangential E-field (V/m) beta 95%-CI p R² Intercept -0.2618 (-0.6676,0.1440) 0.2047 Sex (males) 0.0026 (-0.0082,0.0135) 0.6332 0.00 Age (years) -0.0002 (-0.0009,0.0005) 0.5789 0.00 Montage (ref: C3Fp2) F3P3 -0.0151 (-0.0171, -0.0131) <0.001 0.07 depth (mm) 0.0070 (-0.0048,0.0187) 0.2434 0.01 thickness (mm) 0.0209 (-0.0198,0.0616) 0.3123 0.01 TIV (mm3) 0.0000 (-0.0001,0.0001) 0.5234 0.00 Gyrification 0.0111 (0.0035,0.0187) 0.0043 0.07 17 Fractal Dimension 0.0280 (-0.0755,0.1316) 0.5938 0.00 Random effects Sigma subject 0.0247 0.021193-0.028745 Total R2 0.9617 Fixed effects Normal E-field (V/m) beta 95%-CI p R² Intercept -0.1973 (-0.4408, 0.0462) 0.1116 Sex (males) 0.0003 (-0.0058,0.0073) 0.8206 0.00 Age (years) -0.0001 (-0.0005,0.0004) 0.8644 0.00 Montage (ref: C3Fp2) F3P3 -0.0166 (-0.0184, -0.0147) <0.001 0.20 depth (mm) 0.0039 (-0.0032,0.0109) 0.2808 0.01 thickness (mm) 0.0179 (-0.0066, 0.0423) 0.1507 0.01 TIV (mm3) 0.0000 (-0.0003,0.0002) 0.7700 0.00 Gyrification 0.0081 (-0.0003,0.0002) <0.001 0.08 Fractal Dimension 0.0173 (0.0036,0.0127) 0.5842 0.00 Random effects sigma subject 0.014383 0.01322-0.01692 Total R2 0.9210 3.4.3 Effect on individual stimulation intensities Similarly, an analysis of the effect of anatomical characteristics (independent variables) on supra-threshold stimulation intensities (dependent variable) was carried out through linear mixed-effects models, for which, results are shown in table 6. Our analysis showed that skin, bone, and CSF volumes had a significant impact on stimulation intensity, with larger volumes requiring larger stimulation intensities to reach a threshold of 0.15 V/m at the hand knob region (Supplementary table 8). In line with these results, multiple regression models for each electrode montage attributed to anatomical variables an inter-subject variance of 42% (C3Fp2) and 52% (P3F3) for the supra-threshold anodal intensity (Supplementary table 9). The same analysis was carried out to analyze the effect of cortical surface parameters on the intensities which induced tangential and normal E-field magnitudes larger than 0.15 V/m on the hand knob. Only gyrification showed a tendency towards significance (p < 0.01), with lower gyrification leading to lower stimulation intensities to reach a threshold of 0.15 V/m. This additional information can be found in (Supplementary Table 10-12. Supplementary Fig. 13-14). First-order Sobol indices derived from mixed linear models relating suprathreshold stimulation intensities (normE >= 0.15 V/m) over the hand knob with individual anatomical characteristics are represented in Fig. 6 (see Supplementary table 13 for more details). According to Sobol indices, the most significant source of variance was the electrode montage used (Si = 0.54) followed by skin (Si = 0.08) and bone (Si = 0.07) volumes. Additional Sobol indices for normal and tangential E-fields are shown in (Supplementary Table 14-15 and Figure 15-16). 18 Table 6. Results of linear mixed model analyses for the supra-threshold individual stimulation intensity (mA) as dependent variable. Significant beta coefficients and p-values are shown in bold. Fixed effects Stimulation intensity (mA) [normE >= 0.15 V/m] beta 95%-CI p R² Intercept 0.1695 (-0.1850, 0.6412) 0.3933 Sex (males) -0.0027 (-0.0058, 0.0000) 0.0505 0.01 Age (years) -0.0374 (-0.1008, 0.0137) 0.1743 0.00 Montage (ref: C3Fp2) F3P3 -0.2705 (-0.2382, -0.1696) <0.001 0.33 Skin 0.2186 (0.1498, 0.2474) <0.001 0.16 Bone 0.2417 (0.1397, 0.3170) <0.001 0.09 CSF 0.1230 (0.0751, 0.2092) <0.001 0.04 GM 0.0024 (-0.0829, -0.1226) 0.9611 0.00 WM 0.0560 (-0.0388, 0.1496) 0.1999 0.00 TIV -0.0004 (-0.0009,0.0000) 0.0816 0.00 Random effects Sigma subject 0.6157 0.5417-0.6861 Total R2 0.6742 Figure 6. Bar chart for Sobol indices corresponding to the mixed model for supra-threshold stimulation intensity (normE >=0.15 V/m) with tissues volumes and montage as independent variables 19 4. DISCUSSION In this study, we addressed the effects of two different tDCS electrode montages, C3Fp2 and P3F3, on inter-subject variability for a sample of 98 subjects. Specifically, a modelling workflow has been proposed in order to analyze and compare individual E-field parameters over three different ROIs, the sensorimotor cortex M1-S1, and supplementary motor areas SMA and preSMA. We first compared the E-field spatial distribution across subjects for the two electrode montages analyzed. Second, we compared individual values for E-field magnitude and focality between the two montages, and we generally found significant differences in the three ROIs considered. We subsequently addressed the relationship between E-field parameters and different anatomical characteristics such as head tissue volumes. The main results obtained are discussed below. 4.1 Spatial correlation across subjects at the fsaverage surface The average E-field distributions on the standardized fsaverage cortical surface (Fig. 2a) agree with those expected for the two electrode montages analyzed, as reported in previous modelling studies [29]. Whereas C3Fp2 shows a contralateral E-field distribution which spreads over the left prefrontal and the right orbitofrontal cortex, P3F3 confines the E-field to the left hemisphere, over M1-S1 and prefrontal motor-related regions. The range of inter-subject deviation found in our sample, with maximum ratios of about 25% (Supplementary Fig. 1-2), is similar to that of 20 to 25 % reported in previous modelling works [35], [76]. Regarding their distribution, we found higher SD values in regions stimulated with higher E-fields (Fig. 2a) for the two electrode montages analyzed. Similar inter-subject correlation values were observed in M1-S1 for both montages C3Fp2 and P3F3 (0.66 and 0.68, respectively; Fig. 3b). However, the range of deviation in correlation values (0.19 to 0.83 and from 0.45 to 0.80 for C3Fp2 and P3F3, respectively) illustrates well the high variability in E-field distribution across subjects for both montages. The larger variation found for C3Fp2 (10%) compared with that of P3F3 (5%) is exemplified in Fig. 3a, where a higher number of stripes (i.e., subjects who presented lower correlation values) is observed. Due to their longer distance, the electrodes in C3Fp2 comprise a larger portion of the head, which may lead to a stronger dependency of E-field distribution on individual characteristics. Conversely, the ability of P3F3 montage in confining more the E-field ipsilaterally over a smaller area in the sensorimotor region seems to be more robust against inter-individual differences. Although there is no definite E- field threshold value for neuromodulation effect, some authors have suggested values around 20 0.1-0.2 V/m. The deviation found in our modelling results is within the range 0-0.1 V/m, which may explain the variability in clinical outcomes reported using these montages [77], [78], [79]. 4.2 Comparison of E-field parameters at individual ROIs for C3Fp2 and P3F3 The position and orientation of the electrodes have shown to be fundamental variables in order to manage E-field dosing parameters in terms of E-field magnitude, focality, and direction [80], [81]. According to our results, P3F3 showed significant differences in these parameters with respect to C3Fp2 in the three ROIs analyzed (Fig. 4). Whereas P3F3 led to higher average and maximum E-field values in M1-S1 than C3Fp2 (Fig. 4a-b), the latter induced higher E-field values in supplementary motor SMA (Fig 4e-f) and preSMA (Fig. 4i-j) regions. In line with these results, Salvador et al. reported higher average E-field values in premotor areas for C3Fp2 when compared with other electrode montages [82]. In terms of focality, significant differences between C3Fp2 and P3F3 were neither observed in M1S1 nor in preSMA regions (Fig. 4c,k). In these two areas, a higher variability in focality was observed for C3Fp2. In SMA, the focality was significantly higher for C3Fp2 (Fig. 4g, Table 1), which also showed a lower deviation of 25.1% compared with the 32.7% of P3F3 (Table 1). Therefore, we may conclude that the different E-field distributions between the two montages C3Fp2 and P3F3 in the targeted M1-S1 do not lead to significant differences in focality, at least for the non-focal large-electrodes montages here studied. When analyzing local E-field deviation in the cortical surface, Mikkonen et. al [18] reported higher inter-subject variability for focal montages using high-definition electrodes, which have not been addressed in our study. In addition, subsequent analyses should be conducted to verify whether the divergence in results reported in the literature may also be due to different percentages selected for analyzing focality as well as in the heterogeneity of electrode montages and brain regions analyzed. 4.3 Individual supra-threshold stimulation intensities There exists a high variability in cortical excitability observed experimentally [64], with nonresponder rates reported in the order of 50% for anodal tDCS (C3Fp2) over motor cortex [78], [83]. In our study, defining the RR as the proportion of subjects that elicited an induced E-field larger than 0.15 V/m with stimulation intensities of under 2 mA, commonly used in experimental studies, both electrode montages showed almost 100% of responders for both the mean E-field norm and the tangential component over the hand knob. Higher stimulation intensities were required for C3Fp2 to reach 0.15 V/m in both cases (p < 0.001). However, when assessing the normal component the RR decreased to 20% for C3Fp2, which is more in line with experimental 21 evidence. Noteworthy, stimulation intensities in the order of 2 mA were necessary to reach the neuromodulation threshold defined (En >= 0.15 V/m) for both montages (see Fig. 5). This stimulation intensity, however, is infrequent in experimental studies [64] which may explain the low rate of responders seen to be due to suboptimal stimulation intensity. In addition, intrasubject studies have reported an intensity of 2 mA to be the most effective in generating modulation of cortical excitability in the [0.2,2] mA range [84]. Since the normal E-field component seems to better capture individual experimental variability, at least when considering a threshold of 0.15 V/m, our results may suggest that the neuromodulatory effects seen in practice may be induced by the normal E-field component, responsible for somatic depolarization, and not by the tangential component, which has been suggested to modulate the excitability of interneurons in recent multiscale computational studies [85], [86]. Therefore, normal E-field may be the preferred E-field parameters used to design individualized dosingbased protocols. 4.4 Effects of individual anatomical characteristics on E-field parameters The influence of some anatomical characteristics on the variability of the induced E-fields was analyzed using linear mixed-effects models. The results in Table 4 confirm some outcomes previously described in the literature, such as the inverse association of E-field magnitude on individual CSF [33], bone [87] and skin [37] volumes irrespective of electrode montage. The skin volume was the parameter that most influenced the variability of the E-field magnitude in both montages with partial R2 values of about 0.25. In this sense, our results show the same trend that the one reported by Antonenko et al in [29], who referred a percentage of variance for the general E-field between 75 and 93% for different electrode montages analyzed. In our study, individual tissue volumes in the ROI accounted for 83-88% of the variance in the average E-field. For the maximum E-field this percentage dropped to 63% for C3Fp2 (Supplementary Table 1), which is closer to the result reported by Opitz et al. in [33], where around 50% of the variability was attributed to anatomical factors. One possible source of divergence with these studies is that they used global head volumes, whereas we used skin and bone volumes restricted to the region below the anode, and CSF, GM and WM volumes in the ROI. Additionally, Antonenko et al. found a significant interaction of tissue volumes with age when analyzing subjects between 20 and 79 years. We did not detect significant association of sex and age with the E-field, in contrast to [88], possibly due to the more limited range of age between 32 and 40 years analysed in our study. Decreases in cortical volume have been described as a function of the age of the subjects, being more pronounced in women [89]. Due to the age range of the sample, these 22 effects are not observed in the present study [90] .In terms of focality, a significant association was only found with GM volume across montages (Table 4). When analyzing linear models separate for each electrode montage, we found that GM and CSF volumes could explain 56% of the variability in the E-field for P3F3, while there was no significant association with any anatomical parameter for C3Fp2 montage. This may be in line with our E-field distribution analysis: since P3-F3 confines the E-field to a smaller region, the anatomy should influence more the outcomes. Alternatively, C3-Fp2 spreads more the E-field, probably minimizing the anatomical influence. In agreement with conclusions derived by Antonenko et al. [29], inter-subject variations of the E-field magnitude were highly robust across the two electrode montages in our sample group, with a considerable percentage of the inter-subject variability explained by variations of general anatomical characteristics such as TIV, skin, bone or CSF volume in the ROI. Interestingly, nonfocal montages seem to be more sensitive to such general individual parameters, possibly due to the large brain area covered by the electrodes. In this sense, gyrification, an indicator of the convolution of the cortical surface, showed to have a significant effect on the magnitude of the tangential and normal E-field components (Table 4). Thus, the inter-subject differences in the curvature of the target ROI among subjects leads to a variation in the orientation of the induced E-field on the cortex. This effect can also be seen in the decrease of the mean value of the spatial correlation coefficients of the E-field components, especially for the normal component (Figure 3). This low inter-subject correlation for the E-field components may be associated with the variability of neurophysiological effects according to their reported role in the activation of different elements of the neural architecture [91], [92]. Finally, the electrode montage, as a modulator of the direction of incidence of the E-field, shows a strong association with the magnitude of both components, in line with the results obtained for the E-field norm in M1S1. On the other hand, our results show that skin, bone and CSF volumes are directly associated with individual supra-threshold stimulation intensities, with higher intensities being required for larger tissue volumes, thus supporting previous studies which demonstrated their key role in the design of personalized doses for anodal motor cortex stimulation [81], [93]. Additionally, gyrification seems to be associated with individual stimulation intensity in our analysis. Since local gyrification has been reported to be sensitive to various neuropsychiatric disorders [94], [95], as well as to aging and cognitive function [96], our results suggest that it should be taken into account for experimental and therapeutic designs, especially in study samples with individual of older age and/or suffering from pathologies with significant changes in gyri, such as schizophrenia [97]. 23 4.5 Limitations and future work This study provides a new modelling methodology for the analysis of the E-field in selected ROIs volumes derived from cortical atlases. The main limitation of the study is the lack of experimental validation of the results obtained. We believe that the experimental analysis of the inter-subject variability associated to P3F3 should be further studied in future works addressing both premotor and motor areas. In addition, the limited range of age of the sample of subjects addressed in our study may not account for further anatomical variations related to age, which may considerably increase the inter-subject variability here predicted. On the other hand, the E-field threshold considered (0.15 V/m) has been selected according to previous experimental [98], [99] and modelling studies [59], [100]. However, there is no consensus on a definitive threshold for neuromodulatory effects, which may vary across subjects exhibiting a variety of source dependencies [101]. In addition, the electrical conductivity of tissues is an intrinsic source of variability in computational analysis. SimNIBS default conductivities have been chosen for the tissues involved in our models, despite the fact that inter-subject variability of conductivity values is reported to be present in healthy subjects [102], [103], [104]. This effect may be amplified by considering individuals suffering from pathologies that may alter the dielectric properties of the tissues in addition to the anatomical and morphological characteristics of the cortex. Finally, in order to improve the understanding of the association between the E-field and its components on neurophysiological effects, as well as with experimental results, more comprehensive models detailing the neural architecture and stimulated neuron groups for each montage are needed. 5. CONCLUSION This study analyses the E-field induced on cortical motor network regions of 98 subjects under two different electrode montages, the conventional C3Fp2 and an alternative montage P3F3, using MRI-based computational models of the head and brain. For the three main research questions addressed in our study, our simulated results predicted: 1) A high inter-subject variability in E-fields for both electrode montages, especially for the normal component over C3Fp2. 2) Significant individual differences in E-field magnitude and focality between these two montages in M1S1 and in supplementary motor areas SMA and preSMA. 24 3) Anatomical characteristics such as skin, bone and CSF volumes were significant sources of E-field variability in ROIs, irrespective of electrode montage. Gyrification, an indicator of cortical convolution, was also found to be a source of variability for normal E-fields. 4) Similarly to E-field parameters, stimulation intensities were subjected to variability due to anatomical characteristics. The normal E-field component better captured the individual variability and low rate of responder subjects observed in experimental studies. Further combined modelling and experimental studies are still needed in order to better understand whether the predicted E-field differences across subjects correlate with the variability found in experimental studies, as well as with neurophysiological effects observed. 6. ACKNOWLEDGMENTS This work was supported by the was supported by the Asociación Universitaria Iberoamericana de Postgrado (AUIP), Spanish Ministry of Science and Innovation (RTC2019-007150-1, PID2021- 127034OA-I00), Instituto de Salud Carlos III-Fondo Europeo de Desarrollo Regional (ISCIII-FEDER) (PI16/01575, PI18/01898, PI19/01576, PI20/00613, PI21/01875, PI22/01704), Consejería de Economía, Innovación, Ciencia y Empleo de la Junta de Andalucía (CVI-02526, CTS-7685), Consejería de Salud y Bienestar Social de la Junta de Andalucía (PI-0471-2013, PE-0210-2018, PI- 0459-2018, PE-0186-2019, PIN-0097-2020), Consejería de Transformación Económica, Industria, Conocimiento y Universidades de la Junta de Andalucía (PY20_00896, PY20_00903), OTOVIRT project (PIN-0097-2020) from Andalusian Government co-funded by FEDER and the Fundación Alicia Koplowitz. PMR is member of the European Reference Network for Rare Neurological Diseases (project ID 739510). Acknowledgements to Sofia Rita Fernandes for her exhaustive review of this work. 7. ETHICAL STATEMENT This study has the approval of the ethical committee of the University Hospital Virgen del Rocío (Sevilla, Andalucia, Spain) in whose area the MRIs used in this article were collected. All participants gave written informed consent to participate in the study. The research was conducted in accordance with the principles embodied in the Declaration of Helsinki and in accordance with local statutory requirements. 31 [78] I. Laakso, S. Tanaka, S. Koyama, V. De Santis, and A. Hirata, “Inter-subject variability in electric fields of motor cortical tDCS,” Brain Stimul, vol. 8, no. 5, pp. 906–913, 2015, doi: 10.1016/j.brs.2015.05.002. [79] V. Rawji et al., “tDCS changes in motor excitability are specific to orientation of current flow,” Brain Stimul, vol. 11, no. 2, pp. 289–298, 2018, doi: 10.1016/j.brs.2017.11.001. [80] Á. S. Foerster, Z. Rezaee, W. Paulus, M. A. Nitsche, and A. Dutta, “Effects of cathode location and the size of anode on anodal transcranial direct current stimulation over the leg motor area in healthy humans,” Front Neurosci, vol. 12, no. JUL, pp. 1–12, 2018, doi: 10.3389/fnins.2018.00443. [81] C. Evans, C. Bachmann, J. S. A. Lee, E. Gregoriou, N. Ward, and S. Bestmann, “Dosecontrolled tDCS reduces electric field intensity variability at a cortical target site,” Brain Stimul, vol. 13, no. 1, pp. 125–136, 2020, doi: 10.1016/j.brs.2019.10.004. [82] R. Salvador, C. Wenger, M. A. Nitsche, and P. C. Miranda, “How electrode montage affects transcranial direct current stimulation of the human motor cortex,” Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS, vol. 2015-Novem, pp. 6924–6927, 2015, doi: 10.1109/EMBC.2015.7319985. [83] V. López-Alonso, M. Fernández-del-Olmo, A. Costantini, J. J. Gonzalez-Henriquez, and B. Cheeran, “Intra-individual variability in the response to anodal transcranial direct current stimulation,” Clinical Neurophysiology, vol. 126, no. 12, pp. 2342–2347, 2015, doi: 10.1016/j.clinph.2015.03.022. [84] T. Chew, K. A. Ho, and C. K. Loo, “Inter- and intra-individual variability in response to transcranial direct current stimulation (tDCS) at varying current intensities,” Brain Stimul, vol. 8, no. 6, pp. 1130–1137, 2015, doi: 10.1016/j.brs.2015.07.031. [85] H. Chung, C. Im, H. Seo, and S. C. Jun, “Key factors in the cortical response to transcranial electrical Stimulations—A multi-scale modeling study,” Comput Biol Med, vol. 144, p. 105328, May 2022, doi: 10.1016/J.COMPBIOMED.2022.105328. [86] A. S. Aberra, R. Wang, W. M. Grill, and A. V. Peterchev, “Multi-scale model of axonal and dendritic polarization by transcranial direct current stimulation in realistic head geometry,” Brain Stimul, vol. 16, no. 6, pp. 1776–1791, Nov. 2023, doi: 10.1016/j.brs.2023.11.018. [87] H. Seo, H. I. Kim, and S. C. Jun, “The Effect of a Transcranial Channel as a Skull/Brain Interface in High-Definition Transcranial Direct Current Stimulation - A Computational Study,” Sci Rep, vol. 7, no. December 2016, pp. 1–14, 2017, doi: 10.1038/srep40612. [88] S. Bhattacharjee et al., “Sex difference in tDCS current mediated by changes in cortical anatomy: A study across young, middle and older adults,” Brain Stimul, vol. 15, no. 1, pp. 125–140, 2022, doi: 10.1016/j.brs.2021.11.018. [89] F. Crivello, N. Tzourio-Mazoyer, C. Tzourio, and B. Mazoyer, “Longitudinal assessment of global and regional rate of grey matter atrophy in 1,172 healthy older adults: Modulation by sex and age,” PLoS One, vol. 9, no. 12, Dec. 2014, doi: 10.1371/journal.pone.0114478. [90] J. Seidlitz et al., “Brain charts for the human lifespan,” bioRxiv, pp. 1–34, 2022. 32 [91] S. E. Pépés, A. Draper, G. M. Jackson, and S. R. Jackson, “Effects of age on motor excitability measures from children and adolescents with Tourette syndrome,” Dev Cogn Neurosci, vol. 19, pp. 78–86, 2016, doi: 10.1016/j.dcn.2016.02.005. [92] P. T. Fox et al., “Column-Based Model of Electric Field Excitation of Cerebral Cortex,” Hum Brain Mapp, vol. 22, no. 1, pp. 1–14, 2004, doi: 10.1002/hbm.20006. [93] C. Evans, C. Zich, J. S. A. Lee, N. Ward, and S. Bestmann, “Inter-individual variability in current direction for common tDCS montages,” Neuroimage, vol. 260, no. July, p. 119501, 2022, doi: 10.1016/j.neuroimage.2022.119501. [94] D. Sasabayashi, T. Takahashi, Y. Takayanagi, and M. Suzuki, “Anomalous brain gyrification patterns in major psychiatric disorders: a systematic review and transdiagnostic integration,” Translational Psychiatry, vol. 11, no. 1. Springer Nature, Jun. 01, 2021. doi: 10.1038/s41398-021-01297-8. [95] A. Miola, G. Cattarinussi, M. L. Loré, N. Ghiotto, E. Collantoni, and F. Sambataro, “Brain gyrification in bipolar disorder: a systematic review of neuroimaging studies,” Brain Imaging Behav, doi: 10.1007/s11682-022-00713-x. [96] S. Lamballais, E. J. Vinke, M. W. Vernooij, M. A. Ikram, and R. L. Muetzel, “Cortical gyrification in relation to age and cognition in older adults,” Neuroimage, vol. 212, May 2020, doi: 10.1016/j.neuroimage.2020.116637. [97] D. Sasabayashi et al., “Increased brain gyrification in the schizophrenia spectrum,” 2019, doi: 10.1111/pcn.12939/full. [98] M. A. Nitsche and W. Paulus, “Excitability changes induced in the human motor cortex by weak transcranial direct current stimulation,” Journal of Physiology, vol. 527, no. 3, pp. 633–639, 2000, doi: 10.1111/j.1469-7793.2000.t01-1-00633.x. [99] N. M. Lang A Nitsche W Paulus J C Rothwell R N Lemon, “Effects of transcranial direct current stimulation over the human motor cortex on corticospinal and transcallosal excitability,” Exp Brain Res, vol. 156, pp. 439–443, 2004, doi: 10.1007/s00221-003- 1800-2. [100] P. C. Miranda, A. Mekonnen, R. Salvador, and G. Ruffini, “The electric field in the cortex during transcranial current stimulation,” Neuroimage, vol. 70, pp. 48–58, Apr. 2013, doi: 10.1016/j.neuroimage.2012.12.034. [101] R. Ishii, M. Chauhan, S. Takahashi, T. Sumiyoshi, and Y. Yamada, “Neurobiological Mechanisms of Transcranial Direct Current Stimulation for Psychiatric Disorders; Neurophysiological, Chemical, and Anatomical Considerations,” Frontiers in Human Neuroscience | www.frontiersin.org, vol. 15, p. 631838, 2021, doi: 10.3389/fnhum.2021.631838. [102] H. McCann, G. Pisano, and L. Beltrachini, Variation in Reported Human Head Tissue Electrical Conductivity Values, vol. 32, no. 5. Springer US, 2019. doi: 10.1007/s10548- 019-00710-2. [103] H. McCann and L. Beltrachini, “Does participant’s age impact on tDCS induced fields? Insights from computational simulations,” Biomed Phys Eng Express, vol. 7, no. 4, 2021, doi: 10.1088/2057-1976/ac0547. 33 [104] M. Fernandez-Corazza, S. Turovets, P. Luu, N. Price, C. H. Muravchik, and D. Tucker, “Skull modeling effects in conductivity estimates using parametric electrical impedance tomography,” IEEE Trans Biomed Eng, vol. 65, no. 8, pp. 1785–1797, 2018, doi: 10.1109/TBME.2017.2777143. Dose and inter-subject variability analysis for anodal tDCS stimulation over motor cortical network SUPPLEMENTARY MATERIAL 3. RESULTS 3.1 Comparison of E-field distribution across subjects Supplementary Figure 1. Standard deviation normalized by peak E-field for C3Fp2 montage. Supplementary Figure 2. Standard deviation normalized by peak E-field for P3F3 montage. 3.2 Comparison of individual E-field parameters between montages Supplementary Figure 3. Spatial correlation of individual tangential E-fields over the targeted sensorimotor region across all subjects for two different electrode montages C3FP2 (top left) and P3F3 (top right). Violinplot distribution of the correlation coefficients (bottom). Note that each violinplot includes a boxplot with the median value (white dot), 25% quartile and 75% quartile, upper and lower adjacent values. In addition, outliers are represented as red crosses. Supplementary Figure 4. Spatial correlation of individual normal E-fields over the targeted sensorimotor region across all subjects for two different electrode montages C3FP2 (top left) and P3F3 (top right). Violinplot distribution of the correlation coefficients (bottom). Note that each violinplot includes a boxplot with the median value (white dot), 25% quartile and 75% quartile, upper and lower adjacent values. In addition, outliers are represented as red crosses. normE mean r std median r p75th p25th min max Points outliers CV C3Fp2 0.66 0.08 0.68 0.72 0.62 0.19 0.83 4753 64 11.4 P3F3 0.68 0.04 0.68 0.71 0.65 0.45 0.80 4753 84 6.4 E tangent mean r std median r p75th p25th min max Points outliers CV C3Fp2 0.45 0.06 0.45 0.52 0.40 0.20 0.64 4753 28 13.9 P3F3 0.48 0.05 0.49 0.49 0.45 0.25 0.64 4753 33 11.6 E normal mean r std median r p75th p25th min max Points outliers CV C3Fp2 0.28 0.06 0.28 0.31 0.24 0.06 0.46 4753 35 20.9 P3F3 0.30 0.06 0.30 0.34 0.27 0.11 0.52 4753 34 18.1 Supplementary Table 1. Characteristic parameters of the correlation coefficient distributions for the E-field norm and the tangential and normal components of the E-field. Mean and median, first and third quartiles, standard deviation, coefficient of variation, maximum and minimum values, number of vertices of the region and number of outliers are detailed. 3.3 Effects of anatomical characteristics on individual E-field parameters Supplementary Table 2. Mean (SD) values for anatomical variables such as total intracraneal volume (TIV) and tissue volumes for skin, bone, CSF, GM and WM. TIV (dl) Skin (dl) Bone (dl) CSF (dl) GM (dl) WM (dl) Mean (SD) 14.92 (1.4) 0.26 (0.05) 0.16 (0.02) 0.23 (0.05) 0.37 (0.03) 0.43 (0.05) Supplementary Figure 5. Violinplot distribution of individual anatomical parameters such as skin, bone, CSG, GM and WM volumes (dl). Note that each violinplot includes a boxplot with the median value (white dot), 25% quartile and 75% quartile, upper and lower adjacent values. In addition, outliers are represented as red crosses. Supplementary Figure 6. Simple linear regression models for average E-fields (V/m) and different anatomical parameters such as skin, bone, CSF, GM and WM volume (dl). Beta coefficients and p-values are displayed. Supplementary Figure 7. Simple linear regression models for maximum E-fields (V/m) and different anatomical parameters such as skin, bone, CSF, GM and WM volume (dl). Beta coefficients and p-values are displayed. Supplementary Figure 8. Simple linear regression models for focality (dl) and different anatomical parameters such as skin, bone, CSF, GM and WM volume (dl). Beta coefficients and p-values are displayed. Supplementary Figure 9. Simple linear regression models for average angles (deg) and different anatomical parameters such as skin, bone, CSF, GM and WM volume (dl). Beta coefficients and p-values are displayed. Supplementary Table 3. Results of simple linear regression models for C3Fp2 and P3F3 montages. Dependent variables were: average E-field (V/m), maximum E-field (V/m) and focality (dl) in the M1-S1 ROI. Independent variables were skin, bone, CSF, GM, and WM volumes. Models are adjusted by sex and age. Beta coefficients, 95% CIs, p-values and R2 values are listed. Supplementary Table 4. Results of multiple linear regression models for C3Fp2 and P3F3 montages. Dependent variables were: average E-field (V/m) and maximum E-field (V/m) in the M1-S1 ROI. Independent variables were skin, bone, CSF, GM, and WM volumes. Beta coeficientes, 95% CIs, p-values and R2 values are listed. Dependent variable Montage beta 95% CI p Total R2 Average E-field (V/m) C3Fp2 Intercept Sex_M (males) Age TIV Skin Bone CSF GM WM 0.4010 - - - 0.0004 -0.0380 -0.0541 -0.0320 - - (0.3697, 0.4323) - - - (0.0002, 0.0006) (-0.0437, -0.0323) (-0.0636, -0.0446) (-0.0375, -0.0265) - - <0.001 - - - <0.001 <0.001 <0.001 <0.001 - 0.83 C3Fp2 P3F3 Volume (dl) beta 95% CI p R2 beta 95% CI p R2 Average E- field (V/m) Skin -0.0320 (-0.0432,-0.0208) <0.001 0.40 -0.0403 (-0.0550, -0.0255) <0.001 0.51 Bone -0.0425 (-0.0616,-0.0234) <0.001 0.33 -0.0402 (-0.0588, -0.0217) <0.001 0.47 CSF -0.0258 (-0.0359,-0.0156) <0.001 0.36 -0.0258 (-0.0365, -0.0151) <0.001 0.49 GM 0.0103 (-0.0102, 0.0308) 0.3212 0.20 0.0054 (-0.0161, 0.0269) 0.6185 0.36 WM 0.0086 (-0.0033, 0.0204) 0.1537 0.21 0.0039 (-0.0086, 0.0164) 0.5343 0.37 Max E- field (V/m) Skin -0.0593 (-0.0830,-0.0356) <0.001 0.31 -0.0734 (-0.0990, -0.0480) <0.001 0.53 Bone -0.0762 (-0.1164,-0.0359) <0.001 0.24 -0.0689 (-0.1015, -0.0364) <0.001 0.47 CSF -0.0460 (-0.0675,-0.0245) <0.001 0.27 -0.0258 (-0.0461, -0.0056) 0.0127 0.41 GM 0.0171 (0.0213, 0.8040) 0.4234 0.13 0.0004 (-0.0373, 0.0382) 0.9825 0.37 WM 0.0166 (-0.0078, 0.0410) 0.1797 0.14 0.0075 (-0.0144, 0.0294) 0.4978 0.37 Focality (dl) Skin 0.0183 (-0.0665, 0.1030) 0.6698 0.00 0.10032 (0.0126, 0.1880) 0.0254 0.05 Bone -0.0676 (-0.2097, 0.0746) 0.3480 0.01 0.0667 (-0.0387, 0.1721) 0.2122 0.02 CSF 0.0182 (-0.0570,0.0935) 0.6317 0.00 -0.1197 (-0.1777,-0.0617) <0.001 0.15 GM 0.1582 (0.0477,0.2688) 0.0055 0.08 0.2791 (0.2015, 0.3567) <0.001 0.35 WM 0.0411 (-0.0341,0.1163) 0.0411 0.01 0.1092 (0.0501, 0.1683) <0.001 0.12 P3F3 Intercept Sex_M Age TIV Skin Bone CSF GM WM 0.52268 - - - -0.051823 -0.049204 -0.032227 - - (0.5000, 0.5453) - - (-0.0572, -0.0464) (-0.0583, -0.0401) (-0.0370, -0.0274) - - <0.001 - - <0.001 <0.001 <0.001 - - 0.88 Dependent variable Montage beta 95% CI p Total R2 Maximum E-field (V/m) C3Fp2 Intercept Sex_M Age TIV Skin Bone CSF GM WM 0.7953 - - - -0.0635 -0.0900 -0.0474 - - (0.7282, 0.8624) - - - (-0.0795, - 0.0475) (-0.1169, -0.0628) (-0.0617, -0.0332) <0.001 - - - <0.001 <0.001 <0.001 - - 0.63 P3F3 Intercept Sex_M Age TIV Skin Bone CSF GM WM 0.92985 - - - -0.093619 -0.10398 -0.038655 - - 0.8819 0.9778 - - - -0.1051 -0.0822 -0.1233 -0.0847 -0.0489 -0.0285 - - <0.001 - - - <0.001 <0.001 <0.001 - - 0.83 Supplementary Table 5. Results of multiple linear regression models for C3Fp2 and P3F3 montages. The dependent variable was focality in the M1-S1 ROI. Independent variables were skin, bone, CSF, GM, and WM volumes. Beta coefficients, 95% CIs, p-values and R2 values are listed. Dependent variable Montage beta 95% CI p Total R2 Focality (dl) C3Fp2 Intercept Sex_M Age TIV Skin Bone CSF GM WM 0.6900 - - - - - - - - - (0.6509, 0.7292) - - - - - - - - - <0.001 - - - - - - - - 0.00 P3F3 Intercept Sex_M Age TIV Skin Bone CSF GM WM -0.076885 - - - - - -0.1440 0.3058 - -0.3252 0.1715- - - - - (-0.1863, -0.1017) (0.2412, 0.3704) - 0.032 - - - - <0.001 <0.001 - 0.56 P3F3 Intercept Sex_M Age TIV Gyrification Sulcus depth Fractal dimension Cortical thickness 1.4402 - - - - - - - (1.3892, 1.4912) - - - - - - <0.001 - - - - - - 0.00 Dependent variable Montag e beta 95% CI p Total R2 Intensity (mA) [En >= 0.15 V/m] C3Fp2 Intercept Sex_M Age TIV Gyrification Sulcus depth Fractal dimension Cortical thickness 2.2889 - - - - - - - (2.2172, 2.3607) - - - - - - - <0.001 - - - - - - - 0.00 P3F3 Intercept Sex_M Age TIV Gyrification Sulcus depth Fractal dimension Cortical thickness 1.8800 - - - - - - - (1.8164, 1.9435) - - - - - - <0.001 - - - - - - 0.00 Supplementary Table 12. Results of linear mixed model analyses for morphological parametres of cortical surface: depth (mm), thickness (mm), gyrification (rad) and fractal dimesion. Significant beta coefficients and p-values are shown in bold. Fixed effects Intensity (mA) [Etg >= 0.15 V/m] beta 95%-CI p R² Intercept 2.3248 (0.0915, 4.5581) 0.0414 Sex (males) -0.0443 (-0.1470, 0.0584) 0.3959 0.00 Age (years) 0.0040 (-0.0021, 0.0100) 0.1980 0.00 Montage (ref: C3Fp2) F3P3 0.1105 (-0.0855, 0.1354) <0.001 0.04 depth (mm) -0.0366 (-0.1249, 0.0517) 0.4149 0.01 thickness (mm) 0.0470 (-0.1992, 0.2932) 0.7068 0.00 TIV (mm3) -0.0608 (-0.0002, 0.0006) 0.3250 0.01 Gyrification (rad) -0.0443 (-0.1095, -0.0121) 0.0147 0.05 Fractal Dimension 0.0843 (-0.4452, 0.6139) 0.7537 0.00 Random effects Sigma subject 0.1239 0.0089-0.3046 Total R2 0.9301 Fixed effects Intensity (mA) [En >= 0.15 V/m] beta 95%-CI p R² Intercept 2.7071 (-0.1876, 5.6018) 0.0667 Sex (males) 0.0084 (-0.1247,0.1415) 0.8206 0.00 Age (years) 0.0031 (-0.0047,0.0110) 0.8644 0.00 Montage (ref: C3Fp2) F3P3 0.4090 (0.3596, 0.4584) <0.001 0.26 depth (mm) -0.0001 (-0.1146,0.1144) 0.9983 0.00 thickness (mm) 0.0821 (-0.2369, 0.4012) 0.6122 0.00 TIV (mm3) 0.0002 (-0.0003,0.0007) 0.4220 0.00 Gyrification (rad) -0.0832 (-0.2369,0.4012) <0.010 0.04 Fractal Dimension 0.1657 (-0.5206,0.8520) 0.6345 0.00 Random effects sigma subject 0.3099 0.2716-0.4237 Total R2 0.8613 Supplementary Table 13. Bar chart for Sobol indices corresponding to the mixed model for suprathreshold stimulation intensity (normE >= 0.15 V/m) with tissues volumes and montage as independent variables. Montage Skin Bone CSF GM WM TIV Si 0.54 0.08 0.07 0.02 0.00 0.01 0.00 Supplementary Figure 15. Bar chart for Sobol indices corresponding to the mixed model for suprathreshold stimulation intensity (En >=0.15 V/m) with morphological parameters and montage as independent variables. Supplementary Table 14. Bar chart for Sobol indices corresponding to the mixed model for suprathreshold stimulation intensity (En >= 0.15 V/m) with morphological parameters and montage as independent variables. Montage depth Fractal dimension gyrification thickness Si 0.59 0.00 0.00 0.02 0.01 Supplementary Figure 16. Bar chart for Sobol indices corresponding to the mixed model for suprathreshold stimulation intensity (Et >=0.15 V/m) with morphological parameters and montage as independent variables. Supplementary Table 15. Bar chart for Sobol indices corresponding to the mixed model for suprathreshold stimulation intensity (Et >= 0.15 V/m) with morphological parameters and montage as independent variables. Montage depth Fractal dimension gyrification thickness Si 0.22 0.00 0.00 0.01 0.01