Achieving consistent E-fields across participants through individualized dosing presents a fundamental challenge for tES research. Although participant-specific MRI scans and E-field simulations provide a precise solution, their high cost and limited accessibility hinder widespread use. This paper addressed these constraints by developing robust multiple regression models that could accurately predict peak E-field strength using easily accessible demographic data, head morphology, and inter-electrode parameters. Future experimental validation will be crucial to understanding the impact of our individualized dosing on tES efficacy and treatment outcomes.
Structural MRI scans from the CamCAN repository (available at https://www.mrc-cbu.cam.ac.uk/datasets/camcan/) were used in this analysis. These scans were acquired with a 3T Siemens TIM Trio scanner with a 32-channel head coil at the Medical Research Council Cognition and Brain Science Unit in Cambridge UK. T1w scans were acquired using MPRAGE sequence with the following parameters: TR = 2250 ms, TE = 2.99 ms, TI = 900 ms, FA = 9, field of view (FOV) = 256mm x 240mm x 192mm, voxel size = 1.0mm x 1.0mm x 1.0mm and GRAPPA acceleration factor = 2, TA = 4 min 32 s. T2w scans were acquired with a SPACE sequence with the following parameters: TR = 2800 ms, TE = 408 ms, flip angle = 9, field of view = 256mm x 256mm x 192mm, voxel size = 1.0mm x 1.0mm x 1.0mm, GRAPPA acceleration factor = 2. After processing, the final cohort comprised 418 participants (202 male, 216 female), with a median age of 52 years (range: 18-87).
Electrode positions were selected on the basis of their positional variation on the head and their use in previous works. At low stimulation intensities and frequencies often used in tES, the relationship between the applied current and the resulting E-field strength is linear. Consequently, the E-field generated at a single reference intensity can be directly scaled to predict the E-field strength at other current values. All simulations were therefore based on a reference current of +1mA applied to the anode. In the case of conventional montages with a single cathode this +1mA anode corresponded to a cathode set to 1mA. Conventional montage positions are fronto-central (anode: F4, cathode: Cz), centro-frontal (anode: C3, cathode: FP2), cross-hemisphere frontal (anode: F3, cathode: F4), parietal-frontal (anode: P3, cathode: FP2), cross-hemisphere central (anode: C3, cathode: C4) and midline fronto-occipital (anode: FPz, cathode: Oz). HD montages included one central anode set to 1mA surrounded by four return electrodes set to 0.25mA. The HD montages were positioned as follows: frontal (anode: F4, cathodes: F2, AF4, F6, FC4), central (anode: C3, cathodes: P3, Cz, T7, F3), parietal (anode: P3, cathodes: P1, CP3, P5, PO3) and parietal-occipital (anode: PO8, cathodes: P6, PO4, P10, O10). All montages are visualized in Fig. 6.
E-field simulations were performed using Realistic vOlumetric Approach-based Simulator for Transcranial electric stimulation (ROAST) version 3.0, an open-source, automated, E-field simulation pipeline in MATLAB (2021a). As part of the ROAST pipeline, all structural MRI scans were segmented using SPM12. To produce more accurate meshing than using T1w scans alone, both T1w and T2w MRI scans from all of the 653 participants were segmented into five different tissue types (white matter, grey matter, bone, skin and air). After which, all 653 segmentation maps were manually inspected for imperfections such as scalp contact with the CSF and incomplete eye segmentation. At this stage 200 participants were removed due to low-quality segmentation maps and 254 were manually adjusted to correct for minor errors. 199 segmentation maps were retained without corrections, resulting in a cohort of 453 participants for subsequent analysis.
Following segmentation, electrodes were placed according to the respective montage being simulated. Conventional montages used electrodes with 28mm radii, while HD montages used electrodes with 10mm radii. Both types of electrodes were 1mm thick. Head and electrode meshes were generated with the MATLAB toolbox iso2mesh. Finally, the finite element methods solver, getDP, was used to solve the underlying Laplacian equation. Isotropic electrical conductivities were assigned as follows; white matter S/m, grey matter 0.276 S/m, bone S/m, skin 0.465 S/m, air S/m, gel S/m and electrode S/m. The conductivity value of 0.85 S/m was used for CSF to account for the meninges. In this work, the strength of the E-field was considered, and the E-field direction was discarded.
E-field results were transformed to the MNI152 6th generation reference space. The complete pipeline is illustrated in Figure 7. These transformations used FSL version 6.0 (http://fsl.fmrib.ox.ac.uk/fsl/). Brain extraction (BET), followed by linear transformations (FLIRT) and non-linear transformations (FNIRT), were applied to the T1w MRI scans. 35 participants were excluded at this stage due to the Jacobian determinant of the FNIRT warp fields lying outside of a realistic range of values between 0.01 - 100. The affine transformation matrices and warp fields were subsequently applied to E-field voxels corresponding to the grey matter in the segmentation maps used in the E-field simulation pipeline.
The peak E-field strength, defined as the 95th percentile of the E-field strength within grey matter voxels was extracted. The 95th percentile, rather than the peak voxel, was used to avoid partial volume effects. To ensure data quality, we removed extreme outliers from our dataset. Outliers were defined as peak E-field strength values that fell outside three times the interquartile range (IQR) for each montage. These extreme outliers may have been caused by the non-deterministic nature of mesh generation, which can cause local defects. Of the 4,180 E-fields generated across all 10 montages, 19 were rejected, resulting in the final dataset of 4,161 E-fields.
In addition to E-field strength analysis on focality was conducted on four representative montages. Focality was measured in MNI152 space as grey matter volume with an E-field strength equal to or exceeding 75% of the 99.9th percentile of the E-field distribution, as in previous works. We investigated the relationships between focality and our independent variables using both simple correlational analyses and partial Spearman's rank-order correlations. These findings were largely not significant and are provided in Supplementary Materials Tables 8 and 9.
In addition to age and gender, our tES dosage estimation models used anatomical features such as head circumference, cephalic index, BMI, and inter-electrode distance. The CamCAN repository provided information on age, gender, height, and weight, which were used to calculate BMI. In 54 out of 418 participants, BMI was not calculable due to missing data. For these cases, the missing BMI values were imputed using the median value of 24.05 kg/m (IQR: 20.71-30.56). As head circumferences and cephalic indexes were not included, we estimated these values using MRI scans.
Head circumference was measured by placing 60 points around the transverse plane slice of the MNI152 template; a similar approach was used by Antonenko et al.. The inverse of the participant-specific affine linear transform was then applied to these points, and the sum of the Euclidean distance between the points in native space was used to estimate the head circumference. We compared our measured head circumference values to those predicted by Bushby et al.'s models to validate our approach. This model predicts mean head circumferences of 577mm for males and 553mm for females based on the height of participants in the dataset. This is within 6mm of our approach for males (mean±SD, 583mm±14.5mm) and 1mm of our approach for females (mean±SD, 552 ± 13.5 mm). This suggests that measuring head circumference with a tape measure is comparable to our MRI-based approach.
Head length and width were measured on the same transverse plane slice of the MNI152 template used for head circumference calculation. The length was measured as the distance between the most distal points of the forehead and the back of the head (mean ± SD, 190 mm ± 9.0 mm). The width was measured as the maximum distance between the most lateral points of the head along the axis perpendicular to the length axis (mean ± SD, 158 mm ± 6.7 mm). These measurements were transformed from MNI152 space to native space using the inverse of the linear affine transformation, as with the head circumference calculation. The cephalic index (mean ± SD, 83.1 ± 4.09) is the ratio of the head width to length, calculated as the head width divided by the length and multiplied by 100.
Inter-electrode distances were extracted in each participant's native space after head models were re-oriented to a common right, anterior, superior (RAS) coordinate system. The inter-electrode distance was calculated as the Euclidean distance between the midpoints of the electrode meshes of the stimulation and return electrodes. In conventional montages, the inter-electrode distance was measured between the single return and stimulation electrodes. For HD montages, the inter-electrode distance was the average distance between the stimulation electrode and the four return electrodes.
Finally, data standardization was applied across all data types with appropriate techniques depending on their approximate distributions. Min-max scaling was applied to the uniformly distributed age, and z-score standardization was applied to the normally distributed head circumference, cephalic index, and peak E-field strength. Median interquartile range standardization was applied to the skewed BMI distribution, and gender values were set to 1 (male) or 0 (female). The relationships between these features and the peak E-field strength is shown in Fig. 8.
To train our models to estimate peak E-field strength, we first compiled our testing and training datasets using MRI-informed simulations. We then analyzed the relationship between the features extracted and the peak E-field strengths using simple correlation analyses. Pearson's linear correlation was applied to compare peak E-field strengths against head circumferences and cephalic indices, as these variables were normally distributed. For the non-normally distributed variables, age and BMI, we used Spearman's Rho correlation. Additionally, we performed two-sample t-tests to compare the peak E-field strengths between genders. To understand the independent correlations of each variable, we conducted Spearman's rank-order partial correlations. Bonferroni correction was applied to all tests to adjust for multiple comparisons.
After which, we trained a series of robust MLR models using an iteratively reweighted least squares approach with the bisquare weight function (tuning constant k=4.685). Multicollinearity was assessed across all MLR models using the variance inflation factor where values below 2.5 were considered acceptable. These models were designed to predict peak E-field strength values specific to each montage. The mathematical representation of these montage-specific models is shown in Eq. (1). The predicted montage-specific (MS) peak E-field strength is defined as . weights are present for each predictor, including age (), gender (), head circumference (), cephalic index (), and BMI () as well as the intercept ().
To assess the accuracy of predicted peak E-field strengths () compared to MRI-informed peak E-field strengths when 1mA of stimulation was applied, we employed 5-fold cross-validation. The model's performance on each test set was evaluated using adjusted R and normalized root-mean-square error (NRMSE).
Montage-agnostic models were designed to generalize to unseen electrode configurations, unlike montage-specific models, which are limited to the montage used in their training data. To achieve this, we trained separate montage-agnostic models on HD and conventional montages with 5-fold cross-validation across participants and leave-one-out cross validation across montages.
First, multiple linear regression models were trained for each montage in the training set as before. Next, regression coefficients () were averaged across montages. Then, to incorporate inter-electrode distance which is a differentiating term between montages, E-field magnitudes were scaled () based on the electrode distance term (). This scaling with electrode distance was performed in a linear fashion as shown in Eq. 2.
As well as a non-linear fashion as shown in Eq. 3 to capture the non-linear relationship between electrode distance and E-field strength. As in montage-specific models, model performance on the respective test sets were evaluated using adjusted R and NRMSE.
Given the direct relationship between applied current and resulting E-field strength, the participant-specific current dosage, , can be calculated using Eq. 4. The target peak E-field strength, , can be set to any desired value depending on the specific requirements of a study or application. Traditionally, a current was administered in MRI-informed participant-specific simulations to achieve a generated peak E-field strength, . The values of , , and were then used to determine the participant-specific current dosage, .
In our approach, the proposed models predict when 1 mA is applied (i.e., ). Please note, in the results section, for illustration purposes, we set the target peak E-field strength as 0.1V/m. Levene's tests with Bonferroni correction were applied to compare the variance of the resulting E-field distributions between standardized and fixed dosing.