Personalization of Multi-electrode Setups in tCS/tES: Methods and Advantages

Transcranial current stimulation (tCS or tES) protocols yield results that are highly variable across individuals. Part of this variability results from differences in the electric field (E-field) induced in subjects’ brains during stimulation. The E-field determines how neurons respond to stimulation, and it can be used as a proxy for predicting the concurrent effects of stimulation, like changes in cortical excitability, and, ultimately, its plastic effects. While the use of multichannel systems with small electrodes has provided a more precise tool for delivering tCS, individually variable anatomical parameters like the shape and thickness of tissues affect the E-field distribution for a specific electrode montage. Therefore, using the same montage parameters across subjects does not lead to the homogeneity of E-field amplitude over the desired targets. Here we describe a pipeline that leverages individualized head models combined with montage optimization algorithms to reduce the variability of the E-field distributions over subjects in tCS. We will describe the different steps of the pipeline – namely, MRI segmentation and head model creation, target specification, and montage optimization – and discuss their main advantages and limitations.


Introduction
Transcranial current stimulation (tCS), including transcranial direct current stimulation (tDCS), transcranial alternating current stimulation (tACS), and transcranial random noise stimulation (tRNS), is a family of noninvasive neuromodulatory techniques that employ weak (1-4 mA) electrical currents applied via electrodes placed on the scalp for long durations (20-40 min) [1,2]. Concurrent effects of stimulation range from changes in cortical excitability [3] to modulation of ongoing endogenous oscillations [4]. Hebbian-based mechanisms are hypothesized to lead to long-lasting plastic changes in the brain [5], leading to an increasing interest in putative therapeutic applications in a range of neurological diseases [6]. One factor that limits the usefulness of tCS is the widely reported intersubject variability of responses to stimulation [7]. Several factors can explain variability, but here we will focus on the physical agent of the effects that tCS has on neurons: the electric field (E-field) induced in the tissues.

Biophysical Aspects of tCS
The distribution of currents in the head can be described mathematically by the electric field (E-field) vector induced by tCS. Depending on the orientation of the latter with respect to neuronal processes, the membrane of pyramidal neurons is polarized (approximately 0.2 mV per 1 V/m of E-field value, [8]), which leads to the observed concurrent effects of stimulation. One common hypothesized mechanism is the polarization of the soma of pyramidal cells due to the component of the E-field perpendicular to the cortical surface (E n ), [9]. However, other mechanisms of interaction are possible, such as the polarization of axon terminals [10].
The E-field distribution depends on factors such as head geometry (thickness and shape of the head tissues), electrical properties of the tissues (electrical conductivity, σ), location and geometry of the electrodes, and the currents that are applied via the electrodes [11]. Since in vivo measurements of the E-field still pose a number of technical challenges [12,13] and cannot easily be carried out, computational head models based on structural data (usually head MRIs) are usually employed to estimate it [14,15].
Initial uses of computational head models were limited to a posteriori analysis of the E-field distribution of electrode montages typically applied in experimental protocols [15,16]. In recent years, several algorithms have been described to leverage these head models with the objective of optimizing some dose parameters (position and currents of the stimulating electrodes) to target a specific brain region and/or cortical network [9,[17][18][19][20].
This paradigm shift from "one-model-fits-all" montages to individualized montages leveraging subject-specific head models and dose optimization algorithms has the potential to reduce intersubject variability in the outcomes of tCS and allow for more effective and safe protocols. However, several parameters can affect the outcome of these modeling and optimization pipelines. In this work, we will study how uncertainties in target specification, tissue electrical conductivities, and the threshold for neuromodulatory effects can affect the outcome of the optimization. We will also discuss some of the potential benefits of these pipelines, especially in relation to reduction of intersubject variability of the results of optimization.

Subjects
We included seven healthy children and adolescents (three males) aged 10-17 years (M 14; SD 2). The study was approved by the Ethics Committee of the Faculty of Medicine, Kiel University, Kiel Germany. 1 All participants and their parents were instructed about the study, and written informed consent according to the Declaration of Helsinki on biomedical research involving human subjects was obtained. The study is part of the STIPED project. 2
Tissue segmentation was performed using an in-house implementation combining extra-cerebral tissue segmentations from a new segmentation approach, which will be included in a future version of the open-source simulation toolbox SimNIBS, 3 with brain tissue segmentations and cortical gray matter (GM) surface reconstructions from FreeSurfer [21]. Finite element head models were then generated (see Fig. 1), including representations of the scalp, skull, cerebrospinal fluid (CSF), gray matter, and white matter (WM). The head models also contained representations of Pistim electrodes (1 cm radius, cylindrical Ag/AgCl electrodes 4 ) placed in 61 positions of the 10-10 EEG system. For the electrodes, only the conductive gel underneath the metal connector was represented in the head model. Unless otherwise stated, the scalp, skull, and CSF were modeled as isotropic with conductivities of 0.33 S/m, 0.008 S/m, and 1.79 S/m, respectively, which are appropriate values for the DC-low frequency values used in tCS [15]. The GM and WM were modeled as anisotropic (volume normalization, [22], with isotropic conductivity values used for diffusion tensor scaling of 0.40 S/m-0.15 S/m, for the GM -WM, [15]). E-field calculations were performed in COMSOL 5 using secondorder tetrahedral mesh elements to solve Laplace's equation [11].

Montage Optimization Algorithm
The optimization algorithm used in this study is based on the Stimweaver algorithm [9]. We assumed the normal component of the E-field to the cortical (GM-CSF) surface E n as responsible for the acute effects of stimulation. Positive/negative E n -values, corresponding to E-fields directed into/out-of-the cortical surface, lead to increased/decreased excitability of the soma of pyramidal cells. Inputs to this algorithm include target E n -maps, with information about the target E n -field (E n Target ) in each node of the cortical surface; weight maps, with information about the priority (weight, w) assigned to each node in the optimization; current constraints of the montage (maximum current, in absolute value, per channel I max channel ¼ max i {| I i | } and maximum total injected current I max total ¼ 1 2 P N Channels i¼1 j I i j ); and the maximum number of electrodes in the montage. The objective function in the optimization is the error with respect to no intervention (ERNI, with units of mV 2 /mm 2 ): where N is the number of mesh nodes, N channels is the number of electrodes available for the optimization, E jÀCz n,i is a column vector (lead-field vector) with the normal component of the E-field induced by a bipolar montage that has j as the anode (+1 mA) and Cz as the cathode (À1 mA), and I j is the current (in mA) of electrode j in the montage that is being evaluated. The term P N Channels À1 j¼1 E n,i jÀCz I j yields the normal component of the E-field in the montage being evaluated (for each node i), as follows from the linearity principle [9]. The lead-field terms E jÀCz n,i are calculated on a subject-specific basis using the methods detailed in the previous section. The optimization without the constraint on the number of active electrodes was performed by in-house scripts programmed in Python using the SciPy library. 6 In order to constrain the number of electrodes of the final montage, a genetic algorithm (GA) was implemented following the methods described in [9].

Studies Performed
In this work, we performed several studies to clarify the impact of several inputs and parameters to the optimization algorithm in the results. The first study (study a) aims at determining the impact of target size on the optimization results. This is related to the perceived mechanisms of stimulation underlying the effects of tCS, with some studies focusing on highly localized targets, obtained, for instance, from EEG source localization information [17,20], where other studies focus on more widespread cortical areas with information extracted from cytoarchitectural information [23] or functional imaging data [24]. Optimization algorithms can tackle both of these cases, but it is unclear how current constraints influence the capability of achieving the desired target E n -field with increasing target area size.
Another important parameter that affects the E-field distribution, and therefore the results of the optimization, is the electrical conductivity of the head tissues. Measuring these values in vivo still presents several limitations, and data available in the literature has a wide variability, representing different measuring methods and origin of the tissue samples [25]. Furthermore, some reports indicate that this value might change according to the subject's age, at least for the skull [26]. In study b, we partially tackled this problem by assessing the influence that different conductivity values for the skull tissue have on the optimization results.
In study c, we investigated some of the potential benefits of optimization algorithms on experimental design, namely, less variability on E n -field distribution.
More details about how each study was performed are presented in the next section. In the studies where surface average E n values are presented, they were calculated with the following expression: where A i is the area associated with each node of the mesh (the sum of the areas of all the triangles connected to the node divided by 3) and N patch is the number of nodes in the surface patch where the average is being calculated.

Study A: Effect of Target Size
In this study, we investigated how target size and current constraints affect the actual <E n > achievable on a target given the current constraints. The target was located on the left precentral gyrus, and its size varied from 9 mm 2 (a tiny spot on the gyrus crown) to~275 mm 2 (about the entire left frontal lobe). The different areas were first identified on the cortical surface of a template brain model ( Colin27 7 ), starting from the smallest one and progressively enlarging it up to the maximum size considered. Each area was then remapped separately onto the cortical surface of one of the participants of this study, and single target maps were created: the target area was assigned to excitation, with two values of E n Target (0.25 V/m and 0.50 V/m), and maximum weight w stim ¼ 10; the rest of the cortex was assigned no stimulation (E n Target ¼ 0 V/m) with weight w no-stim varying for each area, in such a way that both conditions have the same relative importance to the ERNI calculation: (w no-stim / w stim ) 2 ¼ Area stim /Area no-stim . Figure 2 shows <E n > on the target as a function of target size. The E n distribution results from optimized montages with an unconstrained number of electrodes, obtained for different combinations of E n Target , maximum current per electrode, and total injected current (I max channel , I max total ) including values exceeding the usual safety limit of 4.0 mA.
For all current constraints and both E n Target values, we observe an overall logarithmic decrease of the <E n > with the target size. This decrease is more rapid for the higher E n Target , and it follows a non-monotonic trend that can be correlated to the cortical curvature of the target area and to the distance between the target area and electrodes. In fact, the significant drop at 108 mm 2 w.r.t. the previous size is likely due to the fact that the ROI is now large enough to comprehend both faces of a sulcus, which have surface normalsand consequently normal electric fieldpointing in opposite directions: once averaged over the whole ROI, this results in a decrease of <E n >. With the next area increase, the ROI extends out of the sulcus, over the two adjacent gyri, and approaches the normal projection of the center of the two closest electrodes, which is reflected in the slow increase of <E n >, reaching a local peak at 656 mm 2 . Further area increases, in the second half of the plot, repeat this pattern, ultimately created by the compounding effect of the gyrification of the target area and the distance from the covering electrodes.
Concerning the current constraints, in Fig. 3, we look separately at the influence of the I max total (3a) and of the I max channel (3b). The shaded area in blue in Fig. 3a represents the value of the maximum <E n > achievable in each target, for both E n Target , normalized with respect to E n Target . This maximum <E n > is obtained with a montage optimization with unrestricted maximum individual and total injected current and is the same for both values of E n Target . As we observe, it also decays logarithmically with the target area, from 95% E n Target for the smallest target to 15% E n Target for the largest. In this case the decay is to attribute utterly to the effect of head anatomy and electrode positions. The figure also shows the <E n > obtained with different current constraints, normalized with respect to the maximum <E n > achievable, for both E n Target . We observe that, for E n Target ¼ 0.50 V/m, only with a I max total ¼ 8.0 mA it is possible to induce in all areas a <E n > at least over 80% of the maximum <E n > achievable. On the other hand, a total injected current I max total of 1.0 mA does not reach even the half of the maximum achievable <E n >, in any target, including the smallest. Moreover, we observe that, as a result of the linearity of I max total and E n Target , for the less stringent condition of E n Target ¼ 0.25 V/m, the exact same relative <E n > on each target area can be achieved with half of the total current. Consequently, I max total ¼ 4.0 mA in this case is sufficient to reach at least 80% of the maximum achievable <E n >. In Fig. 3b, we focus on the current constraints considered in studies b and c. This plot indicates that I max channel modulates the <E n > only up to a given target size (which is smaller for the less demanding condition: E n Target ¼ 0.25 V/m). After this threshold area, the only factor influencing <E n > becomes I max total .

Study B: Tissue Conductivity Values
In this study, we assessed how skull conductivity (σ skull ) values affected the optimization results. We tried three different conductivity values: 0.008 S/m (our standard conductivity value which corresponds to a ratio of scalp-to-skull conductivity of 41), 0.011 S/m (scalp-to-skull conductivity ratio of 30), and 0.041 S/m (ratio of 8). These values cover a wide range of values reported in the literature [26]. For one of the subjects in this study, we calculated the lead-field matrix and performed optimizations with a common target: the left dorsolateral prefrontal cortex (lDLPFC) as identified by Brodmann area 46 [27]. The cortical surface nodes in this area were set to a target  Figure 4a displays <E n > on the lDLPFC as a function of σ skull for the different current and target E n constraints. Average E n values increase nonlinearly with σ skull for every set of constraints except for the less stringent one: target E n of 0.25 V/m with (I max channel , I max total ) ¼ (2.0, 4.0) mA. Figure 4b displays the variation of the total injected current (I total ) in each montage with σ skull . I total tends to decrease nonlinearly with increasing σ skull , except for the most stringent constraint: target E n of 0.50 V/m with (I max channel , I max total ) ¼ (1.0, 2.0) mA. In this case, I total stays almost constant at the highest value allowed (2.0 mA). Figure 5 provides the distribution of E n in the cortical surface for all the optimizations performed in this study. For all optimizations and for the highest σ skull , the position of the electrodes in the optimized montage is very similar across the different sets of constraints, with only the current values being different. For lower σ skull values, and especially in the more stringent optimization constraints, the montages also differ in the electrode positions, often employing bigger separations between the anodes and cathodes.
We also evaluated the change in average E n value that would occur when the optimized montages were evaluated in a model with a different σ skull than the one used to derive the montage (σ skull Eval 6 ¼ σ skull Optim ). This led to changes in average E n values ranging from À45% to 137% of the values obtained when σ skull Eval ¼ σ skull Optim . The average E n values decreased when σ skull Eval < σ skull Optim and they increased otherwise.

Study C: Intersubject Variability
In study c, we investigated the advantages that montage optimization brings in terms of intersubject variability of the E-field distribution. To do this, we performed subject-specific optimizations for six of the subjects in this study. For each subject, the optimization parameters were the same as the ones employed in study b. Each optimization was then evaluated not only on the subject's head model from which it was derived but also in all the remaining head models. For each evaluation, we calculated the average <E n > value on the lDLPFC. To test the homoscedasticity of the different distributions, we used Levene's test. To compare the means of the different groups, we used Welch's t-test. As shown in Fig. 6, the variance of <E n > across subjects was significantly lower when using a personalized montage as opposed to a non-personalized montage (Levene's test, p-value<0.05), except in the most stringent optimization (lowest current constraints with the highest E n Target , Levene's test p-value ¼ 0.73). As is also shown in the figure, no statistically significant difference was found between personalized and non-personalized montages when it comes to the group average of <E n > value across subjects (Welch's t-test p-value>0.651), when the target E n is Fig. 4 Average value of E n (in V/m) in the lDLPFC (a) and total injected current (I max total in mA, b) of the optimized montage as a function of skull conductivity (in S/m). The optimizations were obtained for four different constraints maintained constant. Increasing the target E n leads to a statistically significant higher <E n >, regardless of the current constraints (Welch's t-test p-value<7.6 Â 10 À4 ).

Interplay of Target Size, Cortical Geometry, and Optimization Constraints
When analyzing the influence of target size and optimization constraints, we found an expected decrease of <E n > with the target area. As mentioned before, the non-monotonous nature of the decrease could be attributed to the interplay of different parameters: cortical geometry, positions of the electrodes available for the optimization, and optimization parameters (current constraints and target E nfield). For small targets that do not encompass multiple sulci and/or gyrus, it was possible to achieve even the highest average <E n > value (0.5 V/m) provided enough (total injected) current was set as a limit. Limiting the currents to the safety values used in most studies [28], (I max channel , I max total ) ¼ (2.0, 4.0 mA), even <E n > values of 0.25 V/m cannot be achieved. This depends of course on target position and electrode array. For instance, at the bottom of the sulci under some of the electrode positions available, local maxima have been shown to be created due to the funneling effect of the CSF layer [15]. In these regions, higher <E n > values might be ) are shown next to each group of images (a, b, c, and d). The order of the conductivities of the skull within each group of images is the same (see a). The montages were limited to eight channels. The color scale is common to all plots achievable there with the same current constraints. For these small targets, we also found that the constraint on I max channel can limit the achievable <E n >, with higher values being possible when I max channel is set to the same value as I max total . For larger targets, the maximum achievable <E n > decreases, firstly due to the folded nature of the cortical surface and the electrode distribution (as shown by the rapid decay of <E n > even with unconstrained currents) and then to the current limitations. In particular, I max total is the main limiting factor to <E n >, with the constraint on I max channel not mattering as much. This is expected, as larger targets require the distribution of the current on more electrodes to cover the whole area and achieve the same E n .
Although the effects of having a more dense electrode array available for the optimization were not tested in this study, it is likely that they would be more Fig. 6 Effect of current constraints and target E n value on <E n > calculated on the lDLPFC for personalized (blue) and non-personalized (pink) optimized montages. The top plots (a) show the results for the (I max channel , I max total ) ¼ (1.0, 2.0) mA current constraints, whereas the bottom plots (b) show the results for the (I max channel , I max total ) ¼ (2.0, 4.0) mA constraints. Welch's t-tests were performed for comparisons between the means of the different groups. Levene's tests were performed to test for homoscedasticity between groups with the same constraints beneficial for smaller targets (see also [19]), allowing for higher <E n > to be achieved for the same current constraints.

Influence of Skull Conductivity
The influence of the conductivity of the skull and other tissues on the E-field distribution in tCS is a well-established fact [29,30]. Provided enough current was available to the optimization algorithm, all models reached a similar <E n > value on target. For more stringent constraints (lower currents and/or higher target E n -fields), it might not be possible to maintain a similar <E n > across models (this will again depend on target size). For the latter optimizations, we found that the selected montage employs higher currents and a bigger separation between anodes and cathodes to increase <E n >.
In a more realistic scenario, however, the discrepancy between the subject's skull conductivity and the one employed in the model is the main concern. As our results indicate, this can lead to very big discrepancies between the planned and effective <E n >. These results stress the need for assessing subject-specific tissue conductivity values and use them together with subject-specific computational head models.

Montage Optimization and Intersubject Variability
Consistent with previously published studies [16], we found a large variability when calculating <E n > induced in six head models by non-personalized montages (on average the standard deviation was 21% of the mean value across all cases). Employing personalized montages significantly reduced the variation (standard deviation of 8% of the mean) in all cases except in the more stringent optimization (top-right boxes in Fig. 6). In that case ((I max channel , I max total ) ¼ (1.0, 2.0) mA and E n Target of 0.50 V/m), personalization of montages did not reduce variability or result in an increase in average <E n > at the target. We interpret this as basically showing that the posed optimization problem is very hard to achieve and hence of variable results even with personalization. Increasing the target E n Target to 0.50 V/m does result in a significant increase in <E n > for both current constraints. Ultimately, this may prove to be more important than decreasing the variability of the results. Heterogeneity in <E n > across subjects can always be used as a regressor when analyzing the results of the study (see [31]).

Study Limitations
As all studies involving computational head models, there are a number of limitations in this study related to the simplifications employed by the models. The biggest simplification is the adoption of a homogeneous compartment for the skull tissue, ignoring the spongy bone region [32]. Although this would certainly influence the <E n > values reported here, as well as affect the effects of the different current constraints, it is unlikely that the overall qualitative conclusions of the different studies would be influenced.
Another limitation is related to the fact that we focused exclusively on the normal component of the E-field for optimization. The optimization method employed in this study can easily be used with other E-field components, but the optimized montages would employ electrode positions very different from the ones reported here. Again, this is unlikely to affect the overall conclusions of this study, and its general recommendations can be extended when other components of the E-field are of interest.
Regarding the targets, we only considered connected single target regions, despite the fact that interest has arisen lately regarding applications involving multiple distributed targets (as the ones arising from cortical networks, [5,24]). Again, the stimulation protocol can be readily adapted to these types of targets (with the weights reflecting the statistical significance of the correlation). The conclusions about the influence of current constraints in these types of optimizations are likely to be similar to the ones reported here for the larger areas, but further studies are required.
Finally, we should mention the small number of subjects employed in this study, which limits the generalizability of its results, especially in study c. Future studies are underway which will investigate these findings in a larger population.

Consequences for Protocol Design
The results presented here clearly demonstrate the advantages of employing optimized montages for determining dose parameters in a tCS protocol. They have the potential of reducing variability in the E-field distribution across subjects in a study by taking into effect idiosyncratic subject properties, such as individual head anatomy, electrical properties, and even target location. Of course, these improvements require availability of appropriate data, such as MRI scans with parameters optimized for tissue segmentation [33], protocols for noninvasive determination of tissue electrical conductivities in a fast and reliable way [34,35], as well as a combination of functional and structural data to determine the target for optimization. Regarding target size and location, this should guide the determination of the electrical current constraints of the study, as is clearly illustrated by the previous results.
Another important limitation of the usefulness of montage optimization is the lack of information about the mechanisms of tCS. However, several studies have been published illustrating the interaction of the E-field with neurons [10] and the network amplification effects that can be responsible for the ultimate long-term effects of the intervention [4]. The next step in developing montage optimization protocols would be to combine information about biophysical aspects of current propagation and electrophysiological aspects of E-fieldneuron interaction and neuron-neuron communication [5,36,37].