Individually optimized multi-channel tDCS for targeting somatosensory cortex

OBJECTIVE
Transcranial direct current stimulation (tDCS) is a non-invasive neuro-modulation technique that delivers current through the scalp by a pair of patch electrodes (2-Patch). This study proposes a new multi-channel tDCS (mc-tDCS) optimization method, the distributed constrained maximum intensity (D-CMI) approach. For targeting the P20/N20 somatosensory source at Brodmann area 3b, an integrated combined magnetoencephalography (MEG) and electroencephalography (EEG) source analysis is used with individualized skull conductivity calibrated realistic head modeling.


METHODS
Simulated electric fields (EF) for our new D-CMI method and the already known maximum intensity (MI), alternating direction method of multipliers (ADMM) and 2-Patch methods were produced and compared for the individualized P20/N20 somatosensory target for 10 subjects.


RESULTS
D-CMI and MI showed highest intensities parallel to the P20/N20 target compared to ADMM and 2-Patch, with ADMM achieving highest focality. D-CMI showed a slight reduction in intensity compared to MI while reducing side effects and skin level sensations by current distribution over multiple stimulation electrodes.


CONCLUSION
Individualized D-CMI montages are preferred for our follow up somatosensory experiment to provide a good balance between high current intensities at the target and reduced side effects and skin sensations.


SIGNIFICANCE
An integrated combined MEG and EEG source analysis with D-CMI montages for mc-tDCS stimulation potentially can improve control, reproducibility and reduce sensitivity differences between sham and real stimulations.


Introduction
Transcranial direct current stimulation (tDCS) is a non-invasive brain stimulation method that aims to modulate excitatory or inhi-bitory neural activity in the brain (Antal et al., 2017;Lefaucheur, 2016;Nitsche et al., 2008;Schutter and Wischnewski, 2016). The standard tDCS montage to apply electric currents (mostly 2 mA) on the scalp is by a pair of two large patch-like sponge electrodes (25-35 cm 2 ) (2-Patch). In general, for somatomotor applications (Matsunaga et al., 2004;Nitsche et al., 2003;Nitsche and Paulus, 2000), an anodal patch electrode is placed over the primary motor or somatosensory cortex and a cathodal patch electrode over the supraorbital area, contra-and ipsi-lateral to the side of stimulation, respectively. Due to the broadly distributed electric fields in the brain produced by this so-called anodal stimulation by the 2-Patch montage, tDCS results might suffer from inconsistencies (Veniero et al., 2017), intra- (Antal et al., 2015;Horvath et al., 2015a,b), and inter-subject variability (Laakso et al., 2015;López-Alonso et al., 2014). The cause of the variability might also be attributed to the lack of consideration of an individual targeting and to different conductive profiles of head tissues and anatomical and functional differences among subjects Antonakakis et al., 2019;Huang et al., 2017;Opitz et al., 2015;Wagner et al., 2013). For example, Laakso et al. (2019) found a correlation between the modeled electric field intensity and the efficacy of tDCS in a motor evoked potential experiment, which means that inter-subject variability might be explained by differences in individual electric fields and thus that individual targeting and optimization might improve individual TES efficacy. For an individual targeting, not only target location is relevant, but especially also target orientation. Creutzfeldt and coworkers (Creutzfeldt et al., 1962), who studied the effect of transcortical DC currents in the motor and visual cortex of the cat, showed that neurons are activated by radially-inwards and inhibited by radially-outwards (with regard to the cortical surface) oriented currents. Therefore, anodal stimulation might in fact excite underlying cortical regions, if at least parts of the target area are at radially-oriented gyral crowns or sulcal valleys, while this stimulation might be suboptimal for the mainly tangentiallyoriented targets on sulcal walls (Krieg et al., 2015;Krieg et al., 2013;Radman et al., 2009;Seo et al., 2017). Target areas are also often thought of including excitatory or inhibitory networks, which will thus be parameterized in the terminology of this work by a target orientation that differs by 180°. An appropriate targeting thus means that (1) the injected current should not only be maximal in the target region-of-interest (ROI) in the brain (intensity) and (2) minimal in other areas (focality) but also (3) predominantly oriented parallel (excitation) or anti-parallel (inhibition) to the target orientation (directionality) (Creutzfeldt et al., 1962;Krieg et al., 2015;Krieg et al., 2013;Seo et al., 2017;Wagner et al., 2016). Because of the complexity of such a targeting, multi-channel (mc-) tDCS hardware combined with optimization methods have recently gained considerable interest to achieve an efficient trade-off between intensity, focality, and directionality (Dmochowski et al., 2011;Guler et al., 2016;Ruffini et al., 2014;Sadleir et al., 2012;Wagner et al., 2016;Fernández-Corazza et al., 2020). The mc-tDCS optimization (the tDCS inverse problem) includes the simulation of electric fields in the individual brain resulting from a stimulation at the head surface using a quasistatic approximation of Maxwell's equations (the tDCS forward problem) (Dmochowski et al., 2011;Liu et al., 2018;Miranda et al., 2013;Polanía et al., 2018;Wagner et al., 2013). In this regard, for efficient targeting, the goal is first to determine the target individually and then utilize an appropriate inverse optimization method based on accurate forward simulations, to adapt the mc-tDCS montage individually for each subject, with the goal to achieve an improved neurophysiological stimulation effect in a subsequent tDCS experiment (Dmochowski et al., 2013;Laakso et al., 2019). In this way, differences in target location and orientation among subjects are taken into account. The individualized tDCS inverse approach also needs personalized head volume conductor forward modeling, not only concerning tissue geometries but also to individual tissue conductivities, and the most important conductivity parameter is the one of the skull as found in recent sensitivity investigations (Saturnino et al., 2019;Schmidt et al., 2015).
In this group study with 10 healthy subjects, we will first use combined electroencephalography (EEG) and magnetoencephalog-raphy (MEG) source analysis to reconstruct the main underlying source of the somatosensory evoked potential (SEP) and field (SEF) component at 20 ms post-stimulus, the P20/N20 component. This main source of P20/N20 activity is located in the primary somatosensory cortex (SI) in Brodmann area 3b (Allison et al., 1991;Antonakakis et al., 2020;Antonakakis et al., 2019;Hari et al., 1993;Nakamura et al., 1998). Source analysis will be based on realistic finite element method (FEM) head modeling. Head modeling is personalized not only with respect to the head tissue geometries, but skull conductivity is also estimated individually using an SEF/SEP calibration procedure Antonakakis et al., 2019;Aydin, 2014). This multi-modal approach to reconstruct the P20/N20 component is used to take full advantage of the measured EEG and MEG modalities as they provide complementary information for the same underlying sources. Previous studies have shown in theory (Dassios et al., 2007) and practice (Aydin et al., 2017;Aydin et al., 2015;Aydin, 2014;Fuchs et al., 1998;Huang et al., 2007) that source reconstructions from combined MEG/EEG can outperform single modality one's. A detailed investigation in Antonakakis et al., 2019) has furthermore shown that (i) a combined MEG/EEG approach for the P20/N20 component enables a stable and accurate modeling of not only the source location in Brodmann area 3b, but especially also its orientation, which is not possible when only using a single modality, as well as (ii) taking into account individual skull conductivity variability.
We will also propose a new mc-tDCS optimization method, the distributed constrained maximum intensity (D-CMI) approach, to compute individual stimulation montages for the reconstructed targets. D-CMI includes the concepts of maximum intensity (MI) (Dmochowski et al., 2011) and constrained MI optimization (CMI) (Guler et al., 2016), but it has the additional goal to further distribute the optimization currents and thereby produce less tingling in the skin level. While CMI has already been evaluated in a rehabilitation after stroke study (Dmochowski et al., 2013), D-CMI is presented here for the first time. For specific choices of parameters, D-CMI can be identical to MI or CMI, so that D-CMI unifies and extends the class of intensity-optimization schemes. The proposed new D-CMI mc-tDCS optimization pipeline does not only consider individual targeting (with regard to location and orientation) and head modeling. It also takes into account different experimental parameters such as safety limits, availability of a limited number of stimulation electrodes, limiting the current per electrode, and limiting the skin sensations. Based on MI and CMI, we will present the new D-CMI method and compare it to the alternating direction method of multipliers (ADMM), a mc-tDCS approach for focality-optimization (Wagner et al., 2016), and to the standard 2-Patch method. On the focality-intensity scale (Fernández-Corazza et al., 2020;Dmochowski et al., 2011), the ADMM method (Wagner et al., 2016) used in this study is selected as an approach that represents focality instead of intensity-based tDCS montage optimization. Thus, in contrast to maximum intensity approaches, in ADMM induced currents aim for a focal stimulation of the target area, while minimizing currents in non-target regions. ADMM has already been used in an auditory experiment, where we could show that individualized transcranial electric stimulation increases gap detection performance (Baltus et al., 2018). Our study is preparing a follow-up SEF/mc-tDCS/SEF experiment, where we will then read out the effects of the individualized mc-tDCS stimulation through the differences in pre-and poststimulation SEF amplitudes of early-latency components such as the P20, P22, N30, and P45. Based on the defined quantification metrics, visualization results, and the requirements for this future experiment, we consider our new D-CMI montages to be most appropriate for our goal here of targeting the P20/N20 source in SI using hardware available in our laboratory with only moderate skin sensations. Since in our follow-up experiment, by means of using SEF, we will only read out the stimulation effects on the somatosensory network, we are more interested in intensitybased optimization approaches when considering the intensityfocality scale (Fernández-Corazza et al., 2020).

Subjects
Ten healthy subjects (28 ± 9 years, 8 males and 2 females) participated in this study. The subjects had no history of psychiatric or neurological disorders and had given written informed consent before the experiment.

Somatosensory evoked potential (SEP) and field (SEF) recording and preprocessing
The subject's SEF and SEP were recorded following electrical stimulation of the index finger of the right hand using combined MEG/EEG. The electrical stimuli had a pulse width of 0.2 ms and the inter-stimulus interval varied randomly between 350 ms and 450 ms to avoid habituation. The combined MEG and EEG measurement had 4 runs of 10 minutes each with a sampling rate of 1200 Hz and online lowpass filtering of 300 Hz. The measurement was conducted in a magnetically shielded room (Vacuumschmelze, Hanau, Germany). For EEG measurement, 80 AgCl sintered ring electrodes (EASYCAP GmbH, Hersching, Germany, 74 EEG electrodes plus additional six electrodes to detect eye movements) were used.
Before the measurement, the electrode positions of the EEG cap were digitized using a Polhemus device (FASTRAK, Polhemus Incorporated, Colchester, VT). Polhemus-related electrode digitization errors can influence the accuracy of EEG to MRI registration. A well-known source of error for a Polhemus device is metal in the environment, which has been minimized in our institutional setup (wooden house). We also found that, when using optimal Polhemus reference coil positions close to Cz in combination with our chin rest, registration errors can be kept in the sub-millimeter range, as long as head rotations around this axis are avoided. Such head rotations are difficult to perform for the subjects due to the chin rest. Additionally, we train our technical assistants to specifically take care that such head rotations are avoided. However, the technician's ability to digitize points on the scalp surface repeatedly with accuracy and the variability of measurements between multiple technicians might result in remaining registration errors. In such a case, our postprocessing software allows to project electrodes onto the segmented head surface to correct for remaining Polhemus-related registration errors.
For MEG recording, a whole head system with 275 axial gradiometers and 29 reference sensors (OMEGA 2005, VSM MedTech Ltd. Canada) was used. During the MEG recording, the head position was tracked by three magnetic coils, placed on nasion, left and right preauricular points to determine the subject's head position in relation to the helmet.
After the MEG and EEG acquisition, the SEF/SEP data were preprocessed using CURRY8. 1 Following the pre-processing steps proposed by (Buchner et al., 1994), we applied bandpass filtering between 20 Hz and 250 Hz, notch filtering of 50 Hz (power line noise), and deselected bad channels (EEG) by visual inspection. Trials were cut down to 200 ms duration (50 ms pre-stimulus and 150 ms post-stimulus) and bad trials were rejected using a threshold-based semi-automatic procedure offered in CURRY8 followed by visual inspection of the candidate bad trials in each modality. After this pre-processing, we averaged the remaining approximately 4000 trials to generate the SEF/SEP responses.
The signal-to-noise ratio (SNR) is then calculated based on the method described in (Fuchs et al., 1998). In this method the data is whitened by means of each channel's individual noise level (calculated from the pre-stimulus interval) and resulting in a unit less measure for both MEG and EEG. The transformation in to a common unit is needed for a combined MEG and EEG analysis.

MRI acquisition
3D-T1-weighted (T1w), 3D-T2-weighted (T2w), and diffusionweighted (Dw) MRI datasets were acquired using a MAGNETOM Prisma 3.0 T (Release D13, Siemens Medical Solutions, Erlangen, Germany). T1w scans were conducted with fast gradient-echo pulse sequence (TFE) using water selective excitation to avoid fat shift (TR/TE/FW = 2300/3.51 ms/8°, inversion pre-pulse with TI = 1.1 s, cubic voxels of 1 mm edge length) , T2w scans with a turbo spin echo pulse sequence (TR/TE/FA = 3200/408 ms/90°, cubic voxels, 1 mm edge length) and Dw scans with an echo planar imaging sequence (TR/TE/FA = 9500/79 ms/90, cubic voxels, 1.89 mm edge length), with one volume with diffusion sensitivity b = 0 s/mm 2 (i.e., flat diffusion gradient) and 20 volumes with b = 1,000 s/mm 2 in different directions, equally distributed on a sphere. An additional volume with a flat diffusion gradient, but with reversed spatial encoding gradients was scanned and utilized for susceptibility artifact correction (Ruthotto et al., 2012). During T1w measurement, gadolinium markers were placed at the same nasion, left and right preauricular points for landmark-based registration of MEG/EEG to MRI. All measurements were acquired in supine positioning to reduce head movements and to prevent distorting CSF-brain volume conduction effects due to the brain shift (Rice et al., 2013) that would result from measuring MEG/EEG in a sitting position and MRI in a lying position.
2.4. Source analysis to determine the somatosensory P20/N20 targets In order to perform combined MEG/EEG source analysis (forward and inverse modeling) of the somatosensory P20/N20 targets, in a first step, an individualized skull conductivity calibrated realistic forward model was built for each subject, as described in detail in the following:

Segmentation
T1w and T2w MRIs were used to create a six compartment head model for each subject. The head model consists of the segmented scalp, skull compacta (SC), skull spongiosa (SS), cerebrospinal fluid (CSF), gray matter (GM), and white matter (WM) tissues. In summary, in the first step, the tissues scalp, GM, and WM were first segmented from the T1w MRI. In the second step, the T2w MRI was registered to the T1w scan using an affine registration approach implemented in FSL (Jenkinson et al., 2012), resulting in the T2w_T1w image. Using T2w_T1w, the tissues SC, SS, CSF, and brain were segmented following the steps described in detail in Antonakakis et al., 2019). In the third step, all the segmented tissues resulting from T1w and the registered T2w_T1w were combined to create a head model with six compartments for each subject. To reduce computational complexity, without loss of accuracy for the somatosensory application, the head model was cut in-sufficient distance below the skull following (Lanfer et al., 2012).

Finite element method (FEM) mesh generation
In the next step, a hexahedral FEM mesh with 1 mm mesh size was created from the labeled six compartment model using the freely available software SimBio-VGRID. 2 The voxels from the segmented six compartmented model can be used directly as hexahedral elements. In order to increase conformance to the real geometry and to alleviate the stair-case effects of a regular hexahedral voxel grid, a node shift approach was used to smoothen the compartment interfaces resulting in the final 1 mm geometry adapted hexahedral FEM mesh. Nodes on a two-material interface are moved into the direction of the centroid of the set of incident voxels with minority material, i.e., the material occuring three times or less in the 8 surrounding voxels. The adaptation was calculated using a node-shift of 0.33, ensuring that the interior angles at element vertices were convex and the Jacobian determinant in the FEM computations remained positive. It was shown by (Wolters et al., 2007) that this meshing strategy improves numerical accuracy without increasing computation time and memory usage. With regard to realistic volume conductor modeling, node-shifted instead of regular hexahedra should thus be used for meshing purposes.

Modeling of compartment conductivities
For each subject, we set conductivity values of 430 mS/m for scalp (Ramon, 2004), 1790 mS/m for CSF (Baumann et al., 1997) and 330 mS/m for GM (Aydin, 2014). Skull conductivity was individually calibrated, as described later in Section 2.4.6. For DTI construction and white matter conductivity anisotropy modeling, we performed the following steps: Dw MRI images were first corrected for eddy current and susceptibility artifacts using a reversed gradient approach (Ruthotto et al., 2012). The corrected images were then registered to the T2w image and 3D diffusion tensors were derived following (Jenkinson et al., 2012). In the last step, WM conductivity anisotropy tensors were calculated using an effective medium approach (Tuch et al., 2001) and integrated into the geometry adapted hexahedral FEM model following (Rullmann et al., 2009).

Source space construction
For our dipole scanning inverse source analysis approach, we created a source space in the center of the GM compartment with a resolution of 2 mm and without restrictions to source orientations (no normal-constraint). While the forward problem has a unique solution and high resolution (1 mm, see Section 2.4.2) increases numerical accuracy and alleviates skull leakage effects , we limited the inverse solution to a 2 mm source space resolution. This is due to the limited spatial resolution that can be expected even in combined MEG and EEG source analysis scenarios. We ensured that all sources were located inside GM and sufficiently far away from the neighboring tissue compartments to fulfill the so-called Venant condition, i.e., for each source node, the closest FE node should only belong to elements, which are labeled as GM. It must be fulfilled to avoid numerical problems and unrealistic source modeling for the chosen Venant dipole modeling approach (Vorwerk et al., 2014;Wolters et al. 2007;Medani et al. 2015). Otherwise, monopoles might be induced in the neighboring compartments, i.e., white matter and CSF, which would lead to unrealistic source modeling (Vorwerk et al., 2014).

Leadfield computation
For our final forward modeling solutions, leadfields for both MEG and EEG were calculated using the software SimBio. 3 An isoparametric Lagrangian FEM approach with trilinear basis functions was used. For sufficient computational speed, MEG and EEG leadfield bases (Wolters et al., 2004) and an algebraic multigrid preconditioned conjugate gradient (AMG-CG) solver was used, which has proven to be stable for the considered tissue conductivity inhomogeneities and anisotropies (Lew et al., 2009;Wolters et al., 2002).

Individual skull conductivity calibration and targeting using combined MEG/EEG source analysis
Our combined MEG/EEG reconstruction of the P20/N20 somatosensory targets is interwoven with the estimation of individual skull conductivity. This is important because skull conductivity is considerably varying individually and has a large influence on the EEG (and tDCS) forward problem, while this parameter has nearly no influence on the MEG forward problem Aydin, 2014). Furthermore, overlaid thalamic activity at the P20/N20 peak might still influence EEG, but not MEG topographies and, thereby, reconstructions of the individual Brodmann area 3b thereof (Allison et al., 1991;Antonakakis et al., 2020;Fuchs et al., 1998;Götz et al., 2014;Hari et al., 1993;Kakigi, 1994;Nakamura et al., 1998;Piastra et al., 2020;Rezaei et al., 2020). Therefore, and as proposed in Rezaei et al., 2020;Fuchs et al., 1998;Hari et al., 1993;Kakigi, 1994;Nakamura et al., 1998), the MEG 20 ms peak is used in a first step to determine the P20/N20 target location in SI. After fixing this individual Brodmann area 3b localization, in the next step of the combined MEG/EEG calibration procedure, the EEG P20/N20 peak serves for the estimation of individual skull conductivity in combination with the determination of individual Brodmann area 3b target orientation and strength . Algorithmically, the calibration procedure can be summarized as follows: (a) Define a discrete set of skull conductivities, e.g. C = [c 1 , c 2 , . . .. . ...c n ] (b) For each head model with skull conductivity c i , for i = 1. . .n (b.1) Use a deviation scan at the 20 ms SEF peak to determine the location, a first orientation and first magnitude of the dipole source. (b.2) Keep the location of (b.1) fixed and calculate a second orientation and second magnitude using a least square fit of the 20 ms SEP topography to the fixed source location. (b.3) Keep the location of (b.1) and the orientation of (b.2) fixed and calculate a third magnitude using a least squares fit of the 20 ms SEF topography to the fixed source location and orientation. (b.4) For the calculated dipole of (b.3), calculate the residual variance (RV) to the 20 ms SEP topography. (c) Select the conductivity that gives lowest RV from step (b.4).
In summary, the algorithm uses the complementary information provided by the measured P20/N20 MEG and EEG topographies. In our calibration and targeting procedure, while individually estimating the SC conductivity, to avoid overfitting, for SS conductivity, we used the fixed conductivity ratio of 1:3.6 for SC:SS, following the measurements of (Akhtari et al., 2002).

tDCS forward modeling
For tDCS forward modeling the quasi-static approximation of Maxwell's equations for computing the electric potential is justified, yielding the Laplace equation r Á ðrr/Þ ¼ 0 with r being the conductivity tensor and / the electric potential and inhomogeneous Neumann boundary conditions at the two stimulating electrodes (i.e., À1 mA at a fixed cathode and +1 mA at the anode), and homogeneous one's at the remaining model surface (Miranda et al., 2013;Sadleir et al., 2010;Wagner et al., 2016;Wagner et al., 2014). We use the point electrode model (PEM), yielding sufficient accuracy for the tDCS forward problem regarding practical brain stimulation applications, especially with the small electrodes (PISTIM Ag/AgCl electrodes with a 1 cm radius) we use here (Pursiainen et al., 2018). For the numerical approximation of the Laplace equation, with the exception of the different source and boundary conditions, we use the same FEM head modeling approach as in Section 2.4. From the numerically approximated potential at the nodes, we can then compute the electric field E ¼ Àr/ and the current density J ¼ rEfor each geometry-adapted hexahedral mesh element.

Multi-channel tDCS inverse optimization methods
For appropriate targeting and to make optimal use of the recently developed mc-tDCS hardware, individualized optimization protocols have been developed over the past years, focusing either on optimal intensity or focality, depending on the specific stimulation goal (Dmochowski et al., 2011;Guler et al., 2016;Ruffini et al., 2014;Sadleir et al., 2012;Wagner et al., 2016). To prepare our follow-up somatosensory SEF/mc-tDCS/SEF experiment, our mc-tDCS optimization protocol should also take into account (a) fixed positions of the tDCS electrodes on the scalp (in our case: 39), (b) a fixed maximal number of stimulation electrodes (in our case: 8), as described in Section 2.5, (c) accurate and realistic head volume conductor models and accurate numerical field modeling to solve the tDCS forward problem, and safety aspects such as (d) total injection current, in this study 2 mA, (e) limiting the current per electrode and (f) reducing skin sensations to facilitate sham conditioning. These are all important factors with regard to the available hardware and the necessary safety regulations that have to be fulfilled (Antal et al., 2017;Bikson et al., 2016). The optimization results will be later used in a somatosensory experiment with our Starstim-8 (Neuroelectrics, Barcelona, Spain) mc-tDCS system with a maximum of eight out of 39 possible stimulation electrodes (i.e., 39 holes in the rubber cap into which stimulation electrodes can be inserted). Therefore, we also digitally recorded the m ¼ 1; Á Á Á M with M ¼ 39 possible sensor positions, corresponding to the international 10/20 EEG system, for all subjects with a Polhemus measurement device. The M th electrode is fixed as the reference electrode in all of our forward simulations. This means that it has to carry the sum of all positive or negative currents over all other electrodes so that the overall sum over all M electrodes sums up to zero current (Dmochowski et al., 2011;Guler et al., 2016). The P20/N20 somatosensory source in Brodmann area 3b, reconstructed from combined MEG/EEG data as discussed in section 2.4.6, serves as the individual target for each subject. The general goal is then to find an optimally targeting electrode montage for the individual P20/N20 SI target that is additionally fulfilling the above-described side-constraints. For this purpose, first of all, the superposition principle for a linear combination of all possible current injection patterns from the tDCS can be stated as (Dmochowski et al., 2011). As general rule for notation, uppercase bold letters represent matrices, lowercase bold letters represent vectors and non-bold letters, either upper or lower case, represent scalars.
, the FEM simulated current vector (see Section 2.5) in the j's finite element due to stimulation of the i th electrode pair (i.e., a positive unit current of +1 mA at the i th electrode and a negative unit current of À1 mA at the reference electrode M) and N is the number of hexahedral volume elements in the FE discretization. s 2 R MÀ1 is the applied current vector from the (M-1) non-reference electrodes and e 2 R 3N is the resulting simulated forward modeling solution for the current density, i.e., a vectorvalued quantity e r j À Á 2 R 3Â1 in each finite element. The influence matrix A only has to be computed once by solving (M-1) FEM equation systems, as described in Section 2.5 and implemented in Sim-Bio (Wagner et al., 2014). It can then be used to find the optimal mc-tDCS montage that best fits our stimulation goal, i.e., targeting the individual P20/N20 Brodmann area 3b, as well as fulfilling the additional optimization side-constraints.
Here, we will investigate the mc-tDCS optimization methods alternating direction method of multipliers (ADMM) from (Wagner et al., 2016), maximum intensity (MI) from (Dmochowski et al., 2011), the constrained MI (CMI) (Guler et al., 2016) and, presented for the first time in this paper, the distributed CMI (D-CMI) which incorporates both MI and CMI optimizations. In the following, we will recapitulate ADMM, MI and CMI and formulate the D-CMI method. We will then evaluate the performance of the optimization methods, in comparison to each other and also to a standard 2-Patch montage in order to determine the additional advantage of individualized and optimized mc-tDCS for somatosensory stimulation. On the focality-intensity scale (Fernández-Corazza et al., 2020), ADMM will be our representative approach for optimization of focality (Wagner et al., 2016) (the pink zone on the left side of Fig. 1 in the unification approach of (Fernández-Corazza et al., 2020)), while with MI, CMI and D-CMI, our main focus will be on the representatives of optimization for intensity (the blue zone on the right side of Fig. 1 in the unification approach of (Fernández-Corazza et al., 2020)). The reasoning is that our follow-up somatosensory experiment will only contain short (10 min) mc-tDCS sessions so that a considerable stimulation effect on the somatosensory system, read out by means of a comparison of pre-and post-SEF, might be more probable when using intensity optimization approaches.

Alternating direction method of multipliers (ADMM)
The ADMM method as proposed by (Wagner et al., 2016), on the focality-intensity scale (Fernández-Corazza et al., 2020;Homölle, 2016) is more on the focality than on the intensity side of the scale. It is an optimal control problem for a Laplace equation with Neumann boundary conditions with control and point-wise gradient state constraints. It maximizes the current in the target area and target direction while keeping the current in non-target regions under a given bound. The formulation is given as follows subject to w A target s e A. Khan, M. Antonakakis, N. Vogenauer et al. Clinical Neurophysiology 134 (2022) 9-26 where A target 2 R 3Â MÀ1 ð Þ is the submatrix of A that corresponds to the target area, i.e., if the P20/N20 source was found in element j, and o target 2 R 3 is the orientation of the target source. <.,.> indicates the inner product of the threedimensional vectors. w is a weight allowing high currents in the target region while keeping currents in non-target regions below a threshold e. To ensure convexity of the problem and uniqueness of a minimizer and control the applied currents, an L2 regularization term is introduced to penalize the energy of the applied current and an additional L1 term minimizes the number of active electrodes in the minimization procedure subject to w A target s e with a and b the corresponding regularization parameters. Here we chose the same ADMM parameterization as suggested in (Wagner et al., 2016). We also ensure by rescaling, as also proposed by (Wagner et al., 2016), that the safety constraint with regard to the total injected current (2 mA) is fulfilled. The resulting electric fields are then taken for comparison and analysis. For our goal in this work, ADMM as the representative of the class of focal optimization approaches seems sufficient, but it should be mentioned that first comparisons of ADMM with other focality optimization approaches such as LCMV-beamforming (Dmochowski et al., 2011) and leastsquares or weighted least-squares approaches (Dmochowski et al., 2011) point to the superiority of ADMM with regard to its focality (Homölle, 2016, see his Tables 6.2 and 6.3), surely depending also on the choice of parameters.

Maximum intensity (MI)
Due to different side constraints, on the focality-intensity scale (Fernández-Corazza et al., 2020), the MI method proposed by (Dmochowski et al., 2011) is clearly more on the intensity than on the focality side of the scale. The MI formulation is stated as follows subject to k s k 1 2S Total where s is the current injection pattern with a reference electrode current of (À P MÀ1 m¼1 s m ) that makes sure that the overall injected current always sums up to zero and S Total is the total injected current (2 mA in our case). The maximization of intensity and electroencephalography (EEG) sensors with topographies 20 ms post-stimulation (P20/N20) and realistic head model with somatosensory dipole target from combined MEG/EEG (black cone) for subject S1: (a) MEG sensors and 20 ms SEF topography (b) EEG sensors and 20 ms SEP topography (P20/N20) (c) Six compartment segmented head model with compartments scalp, skull compacta (SC), skull spongiosa (SS), cerebrospinal fluid (CSF), grey matter (GM) and white matter (WM), as also indicated in the legend of the grayscale color scheme, and somatosensory dipole target from combined MEG/EEG (black cone) for subject S1 (d) Segmented head model showing skin surface (light brown), cortical surface (dark and light grey), mc-tDCS cap electrode positions with labels and somatosensory dipole target from combined MEG/EEG (black cone) for subject S1.
in the desired direction at the target is a linear programming problem that can be solved by using the CVX toolbox (Grant and Boyd, 2014).

Distributed constrained maximum intensity (D-CMI)
The D-CMI method is an extension of the constrained maximum intensity (CMI) optimization method which was presented by (Guler et al., 2016). The CMI optimization problem can be stated as follows: subject to k s k 1 2S Total and k s k 1 S maxelec where S maxelec is the maximum current limit per electrode. Since with 2 mA, S Total is kept identical throughout our work, we will refer the CMI optimization approach as CMI (S maxelec ) and as CMI in general. D-CMI, presented here for the first time, aims at achieving high intensity in the target area, similar to MI and CMI, but the optimization function and the side-constraints are chosen in a way that the injected currents are further distributed over multiple electrodes, thus also reducing the sensations at the skin level. In D-CMI we introduce an additional L2 regularization term for the optimization function and the side-constraints are chosen so that both the safety constraint for the total current is fulfilled and, as also proposed for CMI, an upper limit for the current at each electrode is realized by the optimization. The D-CMI optimization problem can then be stated as follows: subject to k s k 1 2S Total and k s k 1 S maxelec where k is an L2 norm regularization parameter that adds strict convexity to the problem with regard to existence and uniqueness of a solution (Wagner et al., 2016) and that can be used to further distribute the current over multiple electrodes. Throughout the paper, we will refer to the D-CMI approach also as D-CMI (S maxelec ; k) to clarify the dependence on its two parameters (as S Total is kept constant at 2 mA throughout this study). Because D-CMI (S maxelec ; k ¼ 0) is identical to CMI (S maxelec ) and D-CMI (S maxelec ¼ S Total ; k ¼ 0) is identical to MI, the D-CMI approach unifies and extends the class of intensity-optimization approaches. In the results section, we will perform a parameter identification study for both S maxelec and k.

Standard 2 -Patch
Additionally, we also compare the mc-tDCS optimization methods with the traditional standard 2-Patch stimulation setup. For this purpose, we simulate for each subject two 5 cm Â 5 cm sponge-like tDCS patches with thickness 4 mm and saline-fluidlike conductivity of 1.4 S/m (Wagner et al., 2014). Following the standard 2-Patch montage as used in (Antal et al., 2017;Matsunaga et al., 2004), for the stimulation of the somatosensory network, the patches were centered at the C3 (anode) and FP2 (cathode) electrode locations which were taken from the digitized Polhemus tDCS cap measurement as explained in Section 2.6. The patches were applied with a total injected current of 2 mA.

Quantification metrics and visualization
The choice of the most appropriate tDCS method for our somatosensory SEF/tDCS/SEF experiment depends crucially on whether sufficiently strong injected currents reach the target area in the direction of the targets to have a maximum effect on the neuronal firing rates (Creutzfeldt et al., 1962;Krieg et al., 2015;Krieg et al., 2013;Wagner et al., 2016). However, while in our future experiment reaching high target intensity is in the foreground, other experimental conditions can easily be imagined, where focality is in the foreground in order to avoid side effects of the stimulation of non-target regions, for example, long time stimulation of epilepsy patients with the goal to reduce seizure frequency and severity (Yang et al., 2020). In addition, we need a method whose parameters give us the flexibility and adaptability to best match the available hardware and the desired results of a tDCS study. Quantification metrics help in this decision process in order to understand the differences between the tDCS methods and their specific advantages and disadvantages. Here, we closely follow the metrics used in (Homölle, 2016;Wagner et al., 2016), namely the averaged current intensity in the target region (IT) where X t is the target grey matter region, |X t | its volume and Â the integration variable. This formulation from (Wagner et al., 2016;Homölle, 2016) is for general extended targets. Here, the target area is only the hexahedral mesh element that contains the reconstructed P20/N20 dipole. For constant A target s over X t , what then remains is IT = |A target s|.
The averaged current intensity in non-target regions (INT) is defined as where X nt is the non-target region (all mesh elements in the brain excluding the target mesh element) and jX nt j its volume. The inner product of the simulated current intensity A target swith the target orientation vector o target indicates the so-called directionality (DIR) Furthermore, we measure the percentage of current intensity that is oriented parallel to the target vector, i.e. parallelity (PAR) and the focality (FOC) of an optimization result We use CURRY8, Matlab 4 (MathWorks, Natick, MA, USA) and SCIRun 5 to visualize the P20/N20 topographies (MEG, EEG) as well as the mc-tDCS electrode montages and the corresponding current densities throughout the brain.

Individualized head modeling
The averaging over trials was used to generate the SEF/SEP responses, resulting in an average SNR of 10 ± 2.93 and 8.07 ± 3.1 for SEF and SEP, respectively, over all subjects.
Exemplarily for subject S1, the P20/N20 SEF and SEP topographies together with the MEG and EEG sensors are shown in Fig. 1  (a) and (b), respectively. Fig. 1(c) presents the six compartment head model segmentation, labeled as scalp, skull compacta (SC), skull spongiosa (SS), cerebrospinal fluid (CSF), grey matter (GM), and white matter (WM). The conductivity anisotropy of WM (Section 2.4.3) and the conductivity of the skull compartment (Section 2.4.6) were modeled individually. For the latter, we used the individual P20/N20 SEF and SEP topographies together with the corresponding six compartment anisotropic head model as input to our skull conductivity calibration procedure (Section 2.4.6), which resulted in individual skull conductivities with a mean and a standard deviation for SC and SS of 7.5 ± 5.4 mS/m and 27 ± 19 mS/m, respectively. Finally, in Fig. 1(d), the 39 Polhemus-measured possible stimulation electrode positions of our Starstim-8 neoprene cap, registered on the surface of the head model, are shown together with their labels and the head model.

Somatosensory targeting
As described in Section 2.4.6, we reconstructed P20/N20 targets for each subject using combined EEG and MEG single dipole scans in the individually calibrated realistic head volume conductor models of Section 4.1. This procedure resulted in an individual dipole target for each subject, localized in Brodmann area 3b in the primary somatosensory cortex SI with predominantly tangential orientation. The P20/N20 somatosensory dipole targets reconstructed from combined MEG/EEG are shown exemplarily (in black) for subject S1 in Fig. 1 (c), (d) and Fig. 2 (c), (d).

Individual parameter identification study for D-CMI
As this study is the first to present D-CMI optimization, an individual parameter identification study was performed for both the regularization parameter k and S maxelec to test their sensitivity to the overall result and identify their best individual choice for later comparison with the competing methods.
D-CMI (S maxelec ¼ 1:5 mA; k): We first fixed S maxelec to 1.5 mA, as this value was found to be overall the most tolerant limit for our subjects without feeling discomfort when using our Starstim-8 system. It should be noted that skin sensations increased in our preliminary experiments when two close by electrodes both carried a 1 mA current. It, therefore, doesn't seem sufficient to only control the maximum current per electrode when trying to minimize skin sensations, even if this parameter is one of the most important, but a better distribution over more electrodes also seems valuable. We then examined a range between 0 and 2000 for the energy penalization parameterk. As the goal was to distribute the currents over the available 8 tDCS electrodes to decrease skin sensations, we selected the k for which the D-CMI (S maxelec ¼ 1:5 mA; k) mainly results in 8 active electrodes already in the first optimization step, as shown in Fig. 2(a), (b), and (d) exemplarily for subject S1. We name the k that produces an 8 electrode montage for each subject as individualized k ind . We show in Fig. 2(a) that, when increasing k (x-axis), the directionality metric DIR (y-axis) for the optimized currents in the P20/N20 target area is quite robust and that, as shown in Fig. 2(b), with increasing k (xaxis), the currents are distributed over more and more electrodes (y-axis). Fork ¼ 0 (no regularization, resulting in 4 active electrodes), the resulting directionality DIR is only 6.04% higher than for k ind ¼ 860(the regularization value that results in 8 active electrodes for subject S1), see Fig. 2(a), while the number of active tDCS electrodes increases from 4 to 8 (Fig. 2(b)). The two optimized montages fork ¼ 0 and k ¼ 860 are visualized together with the head model and the target in Fig. 2(c) and (d), respectively.
For D-CMI (S maxelec ¼ 1:5 mA; k ind ¼ 860) (8 electrodes), the 2 mA total current is spread over three anodes with a maximum current of 0.9 mA injected at electrode CP5 (Fig. 2(d)), while only two anodes are used for D-CMI (S maxelec ¼ 1:5 mA; k ¼ 0) = CMI (S maxelec = 1.5 mA) (4 electrodes), with a maximum of 1.5 mA at electrode CP5, which together leads to a considerable reduction in side effects and related sensations at the skin level such as tingling or pain. This is especially also the case for the cathodes, which are even spreaded over 5 electrodes in D-CMI (S maxelec ¼ 1:5 mA; k ind ¼ 860), three of them more distant, when compared to the only 2 cathodes in CMI (S maxelec = 1.5 mA). Fig. 3 shows the k investigation for all subjects using descriptive statistics with boxplots in (a), (b), and (d) (Campbell, 2021). In Fig. 3(a) we show a relationship between active tDCS electrodes (x-axis) and the necessary k to achieve it (y-axis) and in Fig. 3(b) the relationship between active tDCS electrodes (x-axis) and the resulting DIR (y-axis) for the 10 subjects. While a higher number of active tDCS electrodes requires a higher k (Fig. 3(a)), the DIR measure decreases only minimally as the number of electrodes increases ( Fig. 3(b)). This observation can also be complemented by Fig. 3(c), where the minimal decrease for DIR is shown for an increasingk. In Fig. 3(d) we show the two boxplots for the DIR of D-CMI (S maxelec ¼ 1:5 mA; k ¼ 0) = CMI (S maxelec ¼ 1:5 mA), resulting in 4 active electrodes, and D-CMI (S maxelec ¼ 1:5mA; k ind Þ; where k ind individually varies for each subject, resulting in 8 active electrodes. The boxplot shows that the average difference is only 3%, with a maximum of about 6%. D-CMI (S maxelec ; k ¼ 0) = CMI (S maxelec ): In our second investigation, we fix k ¼ 0 and investigateS maxelec , i.e., the maximum current per electrode, a parameter, which gives us another possibility to increase the number of non-zero stimulation electrodes with decreasing value ofS maxelec . The most interesting is to compute CMIðS maxelec ¼ 0:5), resulting in 8 active stimulation electrodes, and compare it with the 8 electrodes that resulted from D-CMI (S maxelec ¼ 1:5 mA; k ind ) for each subject.
As Fig. 4 shows, the 8 electrode D-CMI (S maxelec ¼ 1:5 mA; k ind ), achieves on average a 10% higher DIR (mean) than the corresponding CMI (S maxelec ¼ 0:5 mA) approach. Importantly, because CMI (S maxelec ¼ 0:5 mA) often leads to multiple closeby electrodes of the maximal 0.5 mA strength, the overall skin level sensations in our experimental tests was quite similar to the competing D-CMI approach, even if the maximal current per electrode for D-CMI (S maxelec ¼ 1:5 mA; k ind ) was often higher.
It should be noted that skin level sensation differences are individual and for the subjects difficult to grasp and that the described differences between D-CMI and CMI with regard to both DIR and skin level sensations are relatively small.
In summary, our individual parameter identification in combination with the preliminary experiments on skin level sensations showed that three arguments are in favor of the D-CMI approach when compared to CMI, namely the slightly higher DIR metric for comparable skin level sensations, the stronger distribution of more distant (with regard to the target) electrodes as well as the additional convexity that is added to the optimization functional by the additional energy penalization term (see also the discussion about the elastic net for ADMM in (Wagner et al., 2016)). Therefore, the D-CMI approach with S maxelec ¼ 1:5 mA and a subject-wise individualized k :¼ k ind value to distribute the total current over the available 8 stimulation electrodes out of the 39 possible ones  in our neoprene cap is used in the following and in our follow-up SEF/mc-tDCS/SEF experiment.

Comparison of the tDCS methods
In our subsequent investigations, we will compare the different proposed optimization methods, first of all between each other, and then also in comparison to the standard 2-Patch approach.
In Fig. 5 we compare the optimized and individualized montages of ADMM (upper row), MI (middle row), and D-CMI (lower row) for three (S1, S2, S3) of the ten subjects, targeted to the individually reconstructed somatosensory P20/N20 SI source (in black) from combined MEG/EEG source analysis. The strength of each anodic and cathodic current is indicated, and additionally colorcoded in red and blue, respectively. ADMM (first row) leads to a rather irregular distribution of anodic and cathodic stimulation electrodes. Main anodic electrodes are over the left posterior (CP5 or C3) and main cathodes over left fronto-central regions (FC1, C1). Main electrodes are often surrounded by electrodes with opposite polarity to improve focality and reduce intensity in nontarget regions. The differences between the subjects are considerable and the maximum current for ADMM in Fig. 5 is 1.54 mA (S3) and, over all subjects, 1.76 mA (not shown here). Due to the L1 norm side-constraint, the MI approach (second row), results in an optimized bipolar montage with only one anode that carries the total injection current of 2 mA, and one corresponding cathode (-2 mA). The distances between anode and cathode are larger than in the ADMM result (obviously for S2 and S3, but also for S1 when considering that the main cathodes in ADMM are between FC1 and C1, while MI only uses FZ). The chosen electrode positions are considerably different between the subjects.
The individualized D-CMI (S maxelec ¼ 1:5 mA; k ind ), now also abbreviated as D-CMI, is comparable to the MI result, using mostly the same main electrodes, with the exception of the main cathode for subject S3. An important difference is that D-CMI currents are distributed over multiple neighbouring electrodes so that the maximum used current is below S maxelec ¼ 1:5 mA, which is in Fig. 5 even only reached for the main anode for subject S3, while the absolute values of all other electrodes are even far below this limit. Especially interesting is that D-CMI in subject S3 distributes over five more distant (with regard to the target) cathodes and reduces the stimulation current from À2 mA in MI to À0.66 mA in D-CMI. ADMM and D-CMI thus use all 8 available Starstim-8 stimulation electrodes and currents with lower amplitude than MI, which only uses 2 of them combined with higher injection currents, resulting in higher sensations at the skin level in our preliminary experiments.
In our next investigation, we will compare the individualized and optimized montages with the non-individualized standard 2-Patch approach. We will also visually analyze the resulting current vector fields in the brain.
In Fig. 6, exemplarily for subject S1, we show the results of the ADMM (first column), MI (second column), D-CMI:=D-CMI (S maxelec ¼ 1:5 mA; k ind ) (third column) and 2-Patch standard approach (fourth column), together with the individualized target (in black). The stimulation montages of all approaches are presented in top ( Fig. 6(a)) and frontal view (Fig. 6(b)). Fig. 6(c) and (d) shows the resulting current vector fields in the brain in a full view over the coronal slice through the target (c) and in a zoomed view at the target side (d). Most importantly, Fig. 6(c) clearly shows that MI and D-CMI reach much higher target intensities than ADMM and 2-Patch, while ADMM outperforms all other approaches with regard to focality, as the intensity in non-target areas is overall much lower. As shown in Fig. 6(d), the individually optimized montages reach high directionality of the injected current vector fields to the target area, while this is not the case for the non-individualized 2-Patch approach. All methods do not produce maximal current vector field amplitudes at the (deeper sulcal) target side, but at more lateral gyral crown areas that are closer to the stimulation electrodes (Fig. 6(c) and (d)).
Finally, Fig. 7 complements Fig. 6 by showing on the x-axes the four different methods and on the y-axes the boxplots from the results of all ten subjects for the metrics IT (a), DIR (b), INT (c), PAR (d) and FOC (e).
The boxplots in Fig. 7 together with Tables 1 and 2, showing mean and standard deviation and statistical analysis results for the examined tDCS methods and quantification metrics, have the goal to strengthen the last statements that could already visually be perceived from the current vector fields for subject S1 in Fig. 6(c) and (d), but now using the defined metrics and in a statistic over all subjects. The effect of the induced electric fields on the dipole target region from the four tDCS methods (ADMM, MI, D-CMI and 2-Patch) was evaluated by employing a one way repeated measures ANOVA (analysis of variance) on the quantification metrics (IT, DIR, INT, PAR and FOC) separately as shown in Table 1 column 6. When necessary the Greenhouse-Geisser correction was used to correct for non-sphericity. Post hoc paired sample t tests were then performed to compare for multiple comparisons between the tDCS methods for each quantification metric separately as showed in Table 2. A P value of less than 0.05 was considered significant for all statistical analyses.
A highly statistical significant effect resulted from ANOVA among the methods ADMM, MI, D-CMI and 2-Patch for the quantification metrics IT (F (3,27) = 18.968, p < .001), DIR (F (3, 27) = 19.028, p < .001), INT (F (3,27) = 39.072, p < .001), and FOC (F (3, 27) = 60.153, p < .001). For the quantification metric PAR (F (3,16) = 6.676, p = .021) only a moderate statistically significant effect resulted. Fig. 7(a) and Table 1 clearly show that the highest target intensity IT is achieved with the MI and D-CMI:=D-CMI (S maxelec ¼ 1:5mA; k ind ) approaches, with only a small advantage for MI (mean of 0.15 A/m 2 for MI versus 0.14 A/m 2 for D-CMI), while 2-Patch and ADMM only achieve means of 0.09 A/m 2 and 0.04 A/m 2 , respectively. From the post hoc test for multiple comparisons between the group means as shown in Table 2 column 2, it is also evident that for metric IT all comparisons showed statistical significant differences between the tDCS methods accept IT D-CMI,MI (P = .375). This also indicates that D-CMI and CMI are similarly performing with regard to the IT metric. Similarly, the highest directionalities DIR are achieved with the MI and D-CMI Fig. 4. Directionality (DIR) boxplots (10 subjects) for the two approaches constrained maximum intensity (CMIðS maxelec ¼ 0:5Þ) with 8 active electrodes and for Distribute constrained maximum intensity (D-CMI (S maxelec ¼ 1:5; k ind )) with 8 active electrodes. The legend showing the characteristics of boxplots (10 subjects, grey dots) as mean (red line), 95% confidence interval (95% CI) (pink) and 1 standard deviation (1 SD) (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) approaches (mean of 0.107 A/m 2 for MI and 0.10 A/m 2 for D-CMI), while 2-Patch and ADMM are at means of only 0.04 A/m 2 and 0.03 A/m 2 , respectively ( Fig. 7(b) and Table 1). Post hoc tests for the methods in DIR, as shown in Table 2 column 3, show statistically significant differences for DIR D-CMI,ADMM (P = .001), DIR D-CMI,2-Patch (P = .002), DIR MI,ADMM (P = .001), DIR MI,2-Patch (P = .003) and non-significant differences for DIR D-CMI,MI (P = .476) and DIR -ADMM,2-Patch (P = .161). The non-significant difference between ADMM and 2-Patch indicate that there is directional similarity between ADMM and 2-Patch. On the other side, the ADMM results in the lowest intensity in non-target regions INT with a mean of only 0.003 A/m 2 , strongly outperforming all other approaches (mean of 0.031 A/m 2 , 0.033 A/m 2 and 0.033 A/m 2 for MI, D-CMI, and 2-Patch, respectively) (Fig. 7(c) and Table 1). It is also shown from the post hoc tests for the INT metric, Table 2 column 4, that there are statistically significant differences for INT D-CMI,ADMM (P < .001), INT MI,ADMM (P = .001), INT ADMM,2-Patch (P < .001) indicating that ADMM is a stimulation method optimized for focality compared to the intensity-optimization methods MI and D-CMI and the standard 2-Patch approach. ADMM therefore also results in the highest focality, FOC, with a mean of 13.2, leaving far behind all other approaches (mean of only 4.66, 3.96, and 2.56 for MI, D-CMI, and 2-Patch, respectively) ( Fig. 7(e) and Table 1). Also evident from the post hoc tests for FOC, Table 2 column 6, are FOC D-CMI, ADMM (P < .001), FOC MI,ADMM (P < .001), FOC ADMM,2-Patch (P < .001). With regard to parallelity (PAR), while the non-individualized 2-Patch approach is only at about 50% with a much larger variability, all individually optimized approaches (ADMM, MI, D-CMI) achieve a mean of about 70%, (Fig. 7(d) and Table 1), i.e., their alignment with the P20/N20 SI target orientations is much better. The post hoc tests for PAR, Table 2 column 5, also show with PAR ADMM,2-Patch (P = .038), PAR D-CMI,2-Patch (P = .019) and PAR MI,2-Patch (P = .023) statistically significant differences between the standard 2-Patch and the three optimization methods. Table 3 shows the highest injected current on an electrode (anode) for each subject for the methods D-CMI and ADMM. The highest injected current for MI and 2-Patch are always 2 mA, the reason why they are not presented in the table. As it can be seen the highest injected current among the subjects for D-CMI is 1.5 mA for subjects S3 and S7 and for ADMM it is 1.764 mA for subject S7. Because of the constraint S maxelec ¼ 1:5 mA for the D-CMI approach, the highest possible injected current cannot exceed 1.5 mA, which is not the case for the ADMM.

Targeted mc-tDCS using MEG/EEG source analysis
In this preparation study for a future somatosensory SEF/mc-tDCS/SEF experiment, we reconstructed the underlying SI sources of the somatosensory P20/N20 components in a group of 10 healthy subjects. We individualized target locations and orienta- Fig. 5. Optimized montages the alternating direction method of multipliers (ADMM), maximum intensity (MI) and distributed constrained maximum intensity (D-CMI) as D-CMI (S maxelec ¼ 1:5 mA; k ind ) for three different subjects (S1, S2 and S3). The montages have been optimized according to the individual location and orientation of the reconstructed P20/N20 target (black dipole). tions using combined somatosensory evoked field (SEF) and potential (SEP) MEG/EEG data in skull conductivity calibrated realistic six compartment head models with integrated WM conductivity anisotropy. While the SEP P20/N20 component might at least in some subjects have an overlaid thalamic potential additionally to the main Brodmann area 3b contribution, the MEG signal at 20 ms post-stimulus is not affected by such too deep and too radial thalamic sources (Allison et al., 1991;Antonakakis et al., 2020;Fuchs et al., 1998;Götz et al., 2014;Hari et al., 1993;Kakigi, 1994;Nakamura et al., 1998;Piastra et al., 2020;Rezaei et al., 2020). Our P20/N20 reconstruction and skull conductivity calibration using combined MEG/EEG, therefore, uses the MEG for the localization of the individual Brodmann area 3b. At the same time, the EEG is then exploited for the estimation of individual source orientation and for skull conductivity calibration. Our experimental goal is thus the stimulation of just a single node of the somatosensory network, which might facilitate the effect evaluations. This is an important difference when comparing the goal of our targeting and optimization procedure with the reciprocitybased optimization of (Dmochowski et al., 2017). Furthermore, in most stimulation studies, targets are usually considered only as location-based targets, i.e., no additional orientation information is used for efficient targeting (Matsunaga et al., 2004). The use of the complementary information from EEG and MEG data in our combined EEG and MEG source analysis together with the individually calibrated skull conductivity Antonakakis et al., 2019) offers the advantage of highlighting the individual differences of the somatosensory P20/N20 SI sources among the subjects. These differences are not only in the target location but also in the target orientation since especially the latter might play an important role and should be taken into account for individual targeting (Creutzfeldt et al., 1962;Krieg et al., 2015;Krieg et al., 2013;Radman et al., 2009;Seo et al., 2017). It should also be mentioned that single modality MEG or EEG reconstructions can lead to considerable differences when compared to combined MEG/EEG for the reconstruction of the 20 ms SEP or SEF component, as shown by (Antonakakis et al., 2019). Our study is motivated by Laakso et al. (2019), who found a correlation between the modeled field intensity and the tDCS efficacy in a motor evoked potential experiment so that individual optimization might also help to better control and especially improve the individual stimulation outcome (Baltus et al., 2018).
We then showed to what extent individually optimized mc-tDCS montages improve targeting with regard to important metrics when compared to the non-individualized standard 2-Patch approach. We used the metrics intensity in target region (IT), directionality (DIR), intensity in non-target regions (INT), parallelity (PAR) and focality (FOC). We expect that an improved performance with regard to these metrics can give us better control in our future somatosensory stimulation experiment. We also modified the maximum intensity (MI) mc-tDCS optimization method (Dmochowski et al., 2011) with control over current per electrode (Guler et al., 2016) by an additional energy penalization term, which we called the distributed constrained maximum intensity (D-CMI) approach, according to our experiment's requirements. These are reduced discomfort such as tingling, pain, itching, and burning sensations (Guarienti et al., 2014;McFadden et al., 2011) and safe stimulation (Antal et al., 2017;Bikson et al., 2016), while keeping nearly highest targeting quality as well as hardwarelimitations (8 active stimulation electrodes out of 39 possible one's in our neoprene cap). Therefore, our study also provides a guideline for the preparation of a controlled mc-tDCS stimulation before its application in an experiment, taking into account the practically most relevant stimulation parameters.  Table 1 Results of a one way repeated measures ANOVA (analysis of variance) showing mean and standard deviation (mean ± SD) for the four methods, ADMM (Alternating direction method of multipliers), MI (Maximum intensity), D-CMI (Distributed constrained maximum intensity) and 2-Patch and their quantification metrics, IT (Intensity in target), DIR (Directionality), INT (Intensity in non-target), PAR (Parallelity) and FOC (Focality). Column 6 shows the statistical results with degrees of freedom (df), f values (F) and P values (P) (*P < 0.05, **P < 0.001).

Metrics
ADMM (mean ± SD) MI (mean ± SD) D-CMI (mean ± SD) 2-Patch (mean ± SD) Statistical effect (df, F, P) A. Khan, M. Antonakakis, N. Vogenauer et al. Clinical Neurophysiology 134 (2022) 9-26 5.2. Comparison of stimulation methods and contribution of D-CMI compared to MI and CMI While maximum intensity (MI), constrained maximum intensity (CMI) and distributed CMI (D-CMI) are, when considering the intensity-focality scale (Dmochowski et al., 2011;Fernández-Corazza et al., 2020;Fernández-Corazza et al., 2015), clearly on the intensity side of the scale with their high IT and DIR metrics, the alternating direction method of multipliers (ADMM) (Wagner et al., 2016) is on the focality side with its high FOC and low INT metrics (Figs. 6 and 7). When compared to the standard 2-Patch approach, since all target optimization approaches ADMM, MI, CMI and D-CMI take both individual target location and orientation into account, they align their injected current vector field much better to the target orientation, resulting in a considerably higher PAR metric (Figs. 6 and 7). The 2-Patch approach is also largely outperformed by MI and D-CMI with regard to the IT, DIR, and even FOC metric, and by ADMM with regard to FOC and INT (Figs. 6 and 7).
For all approaches, the individualization of the targeting and optimization seems important, as shown in Fig. 5 by the large differences in the targets and stimulation montages between subjects and in Figs. 6 and 7 by the much better performance of the individualized optimization approaches.
Our future somatosensory experiment will consist of a prestimulation SEF experiment, followed by an individualized and optimized stimulation of the P20/N20 target in Brodmann area 3b, which is again followed by a post-stimulation SEF experiment. Since we expect that the intensity optimization approaches MI, CMI and D-CMI should be able to generate the largest difference when comparing pre-and post-stimulation SEF, due to their considerably larger DIR metric ( Fig. 7(b)), we will focus only on the contribution of D-CMI when compared to MI and CMI. D-CMI is an extension of MI (Dmochowski et al., 2011) and CMI (Guler et al., 2016) that controls the current per electrode as also proposed by CMI, but with an additional energy penalization term to further distribute the current over multiple electrodes, i.e., with the goal to minimize skin level sensations and high field amplitudes in distant brain areas from the target side. While keeping the important attributes of MI, most importantly a high DIR (mean DIR of D-CMI is only 6.5% smaller than for MI, see Table 1), the D-CMI offers maximal flexibility in controlling and reducing the maximal current per electrode (here, with S maxelec ¼ 1:5 mA, it was reduced by 25% when compared to the MI with 2 mA), complemented by the L2-regularization to further distribute the injected currents over the available stimulation electrodes at especially distant sides (in our case: 8 active electrodes of our Starstim-8 system). While the L2 regularization hardly reduces the DIR (Figs. 2  and 3), it reduced tDCS-induced discomforts such as tingling, pain, itching, and burning sensations (Guarienti et al., 2014;McFadden et al., 2011) in our preliminary experiments. It also adds convexity to the optimization function, which can be an important aspect with regard to the uniqueness of the solution, especially in more complicated targeting situations (see Theorem 3.5 in Wagner et al., 2016). It is shown in Fig. 5 that the D-CMI regularization parameter k takes care that distribution farer away from the target is automatically larger than distribution in the proximity of the target. This is especially visible in the D-CMI result for subject S3 in Fig. 5, where the remaining radial orientation component of the target vector automatically leads to much less distribution at the proximate anodal electrodes side (only 2 anodes) than at the distant cathodes side (6 cathodes). In the extreme case of a fully radial and lateral target, an increase in lambda thus cannot lead to much distribution directly over the target (the proximate closed end of the ''current banana"), while the distribution on the distant side in the contralateral hemisphere will be much bigger (the distant ''pealed end of the current banana"). Due to the regularized and weakened cathodal currents, sensations at the skin level and electric field strength at the brain level in the frontal area will be considerably reduced when compared to MI. In MI, the single À2 mA cathode will lead to significant sensations at the skin level and electric field strength at the brain level also in the frontal area beneath the cathode. In summary, this example shows the power of our new D-CMI approach where only a single additional parameter, the regularization parameter, is used by D-CMI to reduce sensations, while keeping nearly the same high target directionality.
Therefore, we decided to choose the D-CMI (S maxelec ¼ 1:5 mA; k ind ) (8 electrodes) for our future somatosensory SEF/mc-tDCS/SEF experiment. However, it should be noted that skin sensations are difficult to grasp. Therefore, further experiments with a larger group of subjects are needed, for example, a statistical evaluation if it is easier for the subjects to distinguish D-CMI (S maxelec ; k) or D-CMI (S maxelec ; k ¼ 0) = CMI (S maxelec ) from sham. Our implementation is flexible; it also allows a parameter adaptation, for example, when anesthetizing the area under the stimulation electrodes to further reduce tDCS-induced discomforts such as tingling, pain, itching, and burning sensations (Guarienti et al., 2014;McFadden et al., 2011).

Potential of D-CMI in individualized multi-channel transcranial electric stimulation (mc-TES)
We believe that the D-CMI method has the potential to improve the effects of transcranial electric stimulation (TES) (Antal et al., 2017) in general, including, besides tDCS, also transcranial alternating current stimulation (tACS) (Baltus et al., 2018;Kasten et al., 2019), transcranial random noise stimulation (tRNS) (Looi et al., 2017;Splittgerber et al., 2020) and transcutaneous spinal direct current stimulation (tSDCS) (Kuck et al., 2017). D-CMI could help to better control important experimental parameters and thereby also contribute to better reproducibility of TES results. The high DIR values in the target areas with the D-CMI method should also result in much better focality when only considering normal-to-cortex components instead of the modulus, as presented in Fig. 7(c, d).
From our results, we can conclude that individually-targeted multi-channel optimized montages together with individualized head modeling should be incorporated in brain stimulation research in order to increase the chance of achieving clearer and more consistent neurophysiological effects, in agreement with recent literature on that topic (Dmochowski et al., 2013;Fischer et al., 2017;Laakso et al., 2019). As shown for example in Fischer et al. (2017), stimulation with a multi-channel montage increased the M1 excitability compared to a classical tDCS montage when targeting a single brain ROI. Interestingly, inter-subject variability still persisted in all the mc-tDCS optimization approaches also including D-CMI, even if at least the variation (standard deviation) in parallelity was considerably reduced by ADMM, MI and D-CMI when compared to the standard 2-Patch approach (Fig. 7(d)). CMI also showed the same inter-subject variability in Fig. 4. We expect that a part of this variability might be alleviated by the use of a denser electrode array (Guler et al., 2016), especially an increase in the number of possible electrode positions within the optimization process, whereas we were limited here by the only 39 openings of our neoprene cap. Another part of the variability might be ''real", due to for example a deeper target, a lower skull conductivity, or a thicker CSF compartment of a particular subject. In this case, our pipeline has the potential to predict and thus interpret the inter-subject variability of stimulation effects in the later brain stimulation experiment, as also recently proposed in (Kasten et al., 2019).
An important additional advantage of the D-CMI is to ease the use of an experimental sham condition. By limiting the maximal current per electrode and further distributing injected currents over multiple electrodes especially at distant sides, the D-CMI will reduce the number of uncomfortable sensations that can occur beneath the electrodes during stimulation such as tingling, pain, itching, and burning sensations, and thereby reduces the sensation difference to the sham condition. Therefore, the use of D-CMI, or MI or CMI combined with a local anesthetics under the electrodes (Guarienti et al., 2014;McFadden et al., 2011), will facilitate the setup of experiments that involve a sham condition. In this way, possibly complemented by an ''Active-Sham" condition as proposed in (Neri et al., 2020), a controlled and consistent sensation throughout the experiment could be achieved.

Limitations
Due to the maximum principle (Wagner et al., 2016), none of the presented stimulation approaches is able to generate a peak intensity at a deeper target side, intensity maxima are always at the closest cortical areas to the stimulation electrodes ( Fig. 6(c), (d)). For non-invasive deep brain stimulation, other technologies are therefore needed such as temporally interfering electric fields (Grossman et al., 2017).
When considering the intensity-focality scale (Dmochowski et al., 2011;Fernández-Corazza et al., 2020;Guler et al., 2016;Homölle, 2016), in this work we compared D-CMI approaches from the intensity side of the scale with ADMM, one of the focalityoptimizing approaches, as well as a standard 2-Patch approach. We only used ADMM as the representative on the focality side of the scale, since intensity optimization with minimized skin sensations was here in the foreground. Furthermore, first comparisons of ADMM with other focality optimization approaches such as LCMVbeamforming (Dmochowski et al., 2011) and least-squares or weighted least-squares approaches (Dmochowski et al., 2011) pointed to the superiority of ADMM with regard to focality (FOC) (Homölle, 2016, see Tables 6.2 and 6.3). A comparison with even more optimization methods, for example (Dmochowski et al., 2011;Guler et al., 2016;Ruffini et al., 2014;Sadleir et al., 2012;Saturnino et al., 2019), would have gone beyond the scope of this work, as the main goal of our work here is the preparation of individualized targeted (constrained) maximum intensity optimized mc-tDCS montages for the group of subjects later also participating in the follow-up SEF/mc-tDCS/SEF experiment.
Until now, for comparison of skin level sensation differences between D-CMI and CMI approaches, we only ran a first test with one co-author of this work, which pointed us to a slight superiority of the proposed D-CMI when compared to CMI, especially with regard to the frontal electrodes. For a necessary statistical proof, an own ethics proposal and specific experiment in a larger group of subjects is needed that was out-of-scope here. However, it is obvious that an additional energy penalization reduces skin level sensations and thereby brain level electric fields, while our simulations here showed that the also obvious reduction of DIR seems negligible.
With regard to forward modeling, our individualized head models for source analysis and mc-tDCS optimization used six tissue compartments, and the conductivity of the skull was individually estimated, as skull conductivity is the most influential conductivity parameter for both EEG source analysis (Vorwerk et al., 2019) and TES simulation (Saturnino et al., 2019;Schmidt et al., 2015) and is considerably varying inter-individually . We can, however, not exclude that possible inter-subject variability for example in skin conductivity will also influence our presented results (Saturnino et al., 2019;Schmidt et al., 2015;Vorwerk et al., 2019). However, due to the insensitivity of MEG to skin conductivity, at least our source localizations should mainly not be affected. The brain skull interface does not only contain CSF (Jiang et al., 2020), but also the meninges (dura matter (Ramon et al., 2014), arachnoid mater, and pia mater) as well as blood vessels (Fiederer et al., 2016). Therefore, even if first simulations show that our SEF/SEP skull conductivity calibration procedure can compensate at least for parts of these individual modeling inaccuracies, the accuracy of the forward modeling should be further improved.
While this study focused on the tDCS methods comparison and efficient targeting by using combined measurement modalities of MEG and EEG in a group of subjects to prepare the follow-up SEF/mc-tDCS/SEF experiment, possible effects of the electrodeelectrolyte spatial mismatch (Chen et al., 2019) and electrode displacement (Ramaraju et al., 2018) were ignored. However, these important aspects should be considered in future simulations and comparisons of different mc-tDCS optimization methods. On the other side, even when regarding that the state of the brain can considerably affect tDCS outcomes (Li et al., 2015), we think that it might be less important for our planned experiment, because early somatosensory SI components can be considered to be exogenous (Buchner et al., 1994;Riitta Hari and Aina Puce, 2017;Matsunaga et al., 2004).
Finally, factors such as electrode shape and size, electrode-skin contact impedance and electrode shunting effects can also influence field distributions in the brain (Pursiainen et al., 2017;Saturnino et al., 2015). However, these effects were investigated using a complete electrode model (CEM) approach both in EEG source analysis (Pursiainen et al., 2017;Saturnino et al., 2015) as well as in tDCS (Laakso et al., 2017;Pursiainen et al., 2018) and found that for the practical application of EEG and tDCS and for our electrode size, the point electrode model (PEM) used here might be sufficiently accurate.

Conclusion and future work
In this work, we motivated the use of combined MEG/EEG source analysis followed by D-CMI in calibrated realistic head models to target the Brodmann area 3b sources of the somatosensory P20/N20 components in a group of ten healthy volunteers. Our work is a preparatory simulation study before the follow-up SEF/mc-tDCS/SEF experiment, taking into account the most important stimulation parameters such as high target directionality (DIR) with lower skin sensations and electric field amplitude in distant brain areas. We used the complementary information from MEG and EEG for an accurate targeting of Brodmann area 3b with regard to both location and especially also orientation. Concerning the ongoing debate between focality and intensity, we addressed this issue by comparing ADMM (high focality) with MI (high intensity two electrode montage) and D-CMI (high intensity multi-channel montage) methods. A further comparison with a standard nonindividualized 2-Patch approach showed the potential of the individual targeting and optimization combined with state-of-the-art multi-channel tDCS hardware. The fulfilment of comfort and safety aspects and the performance with regard to important metrics showed the potential of the D-CMI method for more precise and more consistent experimental outcomes and for facilitating experimental sham-conditioning.
In a future mc-tDCS experiment, we will apply the calculated D-CMI montages and compare their stimulation effects on the somatosensory network with a standard 2-Patch condition for targeting the P20/N20 somatosensory Brodmann area 3b. Our aim is thus to quantify and compare the stimulation effects and test our hypothesis of a clearer and more consistent experimental outcome of the individualized D-CMI approach. This two-step approach of simulation and experimental work we believe has the potential to significantly improve the understanding of the neurophysiological mechanisms involved in tDCS over the somatosensory cortex.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.