Simulating pad-electrodes with high-definition arrays in transcranial electric stimulation

Objective. Research studies on transcranial electric stimulation, including direct current, often use a computational model to provide guidance on the placing of sponge-electrode pads. However, the expertise and computational resources needed for finite element modeling (FEM) make modeling impractical in a clinical setting. Our objective is to make the exploration of different electrode configurations accessible to practitioners. We provide an efficient tool to estimate current distributions for arbitrary pad configurations while obviating the need for complex simulation software. Approach. To efficiently estimate current distributions for arbitrary pad configurations we propose to simulate pads with an array of high-definition (HD) electrodes and use an efficient linear superposition to then quickly evaluate different electrode configurations. Main results. Numerical results on ten different pad configurations on a normal individual show that electric field intensity simulated with the sampled array deviates from the solutions with pads by only 5% and the locations of peak magnitude fields have a 94% overlap when using a dense array of 336 electrodes. Significance. Computationally intensive FEM modeling of the HD array needs to be performed only once, perhaps on a set of standard heads that can be made available to multiple users. The present results confirm that by using these models one can now quickly and accurately explore and select pad-electrode montages to match a particular clinical need.


Introduction
Transcranial direct current stimulation (tDCS) is a noninvasive neurotechnology, which applies small constant currents (ranging from 0.2 to 2 mA) to the surface of the scalp to achieve neural modulation [1,2]. tDCS has shown promise in treatment for neurological disorders such as epilepsy [3], depression [4], Alzheimer's disease [5], Parkinson's disease Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. [6], pain [7] and stroke [8][9][10]. It has also been shown to improve cognitive functions, such as memory and learning in healthy individuals [11]. There has also been a recent increase in research on alternating current stimulation (tACS) to affect cognitive function [12,13].
To estimate current distributions in the brain and to select a proper pad configuration, high-resolution computational models of the human head are often used in clinical research. Current flow models for tDCS are equally applicable for ac stimulation or random noise stimulation so we refer to these collectively as tES. Intricate folding of the cortical surface, varying skull thickness and distribution of cerebrospinal fluid (CSF) can have a large effect on the current flow inside the head [14,15]. To include these details it is important to use anatomical head models of at least 1 mm resolution. The necessary spatial resolution for modeling requires complex specialized software and substantial computational resources [16,17], making accurate finite element modeling (FEM) modeling often impractical. This is a problem regardless of the specific modeling choice (anisotropy, number of tissues modeled, etc). Fast modeling solutions that have been proposed to-date do not include CSF or detailed cortical folding [18,19]. Thus, solutions with the required spatial resolution are not practical for routine use [20].
In this paper we propose to simulate pad-electrodes (5 × 5 cm 2 , 5 × 7 cm 2 ) with a subset of high-definition (HD) electrodes (4 mm radius) that can be carried out in advance. With this approach a given pad can be quickly simulated as the sum of preexisting solutions computed with high spatial resolution. For each HD electrode in the array simulating the pad a FEM is precomputed one time only. The simulation of the pad is then based on an efficient linear superposition of these individual solutions, which can be readily evaluated for different configurations. The present objective is to make the process of exploring different electrode configurations more efficient by separating it from the computationally demanding FEM modeling.
We envision the FEM modeling of the HD array to be performed once on a set of standard heads, while a large number of clinical researchers and practitioners can then use these solutions to explore and select pad-electrode montages depending of their specific clinical target. In the future, when automated segmentation and modeling becomes more routine these precomputed solutions could also be generated for individual patients so that practitioners may adjust electrode placement for individual subjects.
Here we compare the proposed method with models of continuous pad-electrodes and report deviations in terms of field magnitude and spatial distribution for a set of conventional electrode montages on a reference head that has been previously published [21,22].

Model preparation
A head model with different tissue types was constructed from data by a T1-weighted MRI scan from a healthy subject (male, 36 years old). First, automated segmentation was performed using SPM8 (SPM8, Wellcome Trust Centre for Neuroimaging, London, UK) to demarcate the MRI image into six categories: skin, skull bone, CSF, gray matter (GM), white matter (WM), and air [17]. An in-house MATLAB (Mathworks, Natick, MA, USA) script was used to correct the automatic segmentation errors [17]. Residual segmentation errors were finally fixed manually in ScanIP (Simpleware, Exeter, UK). Small HD electrodes (4 mm radius to prevent overlap) were placed automatically on 336 locations on the scalp using a MATLAB script [17]. The montage follows a conventional concentric system with additional electrodes. Location coordinates for 258 electrodes were used following the concentric arrangement used by BioSemi (BioSemi B.V., Amsterdam, Netherlands). Additional rows of electrodes were placed to potentially allow stimulation of deeper or lowerlying cortical targets. We also placed four additional electrodes around the neck which may serve to emulate distant reference electrodes. In total this adds up to 336 electrodes (258 in original concentric system + 74 in additional rows + 4 neck electrodes) (compare figure 1(a)). After generating the six tissue masks and two masks for electrodes and gel, a complete model of the head can be created.

FEM model generation and calculation
To model current flow we assume the head as a volume conductor comprised of distinct tissues with a specific and uniform conductivity (assuming isotropic conductivity for simplicity). Current flow is governed by Laplace's Equation, which has an unique solution if the electric field is continuous at tissue boundaries and specified at the outer boundary of the volume. For arbitrary shaped media, no closed form solution exists and one typically uses numerical techniques. Therefore, the volume is discretized into a set of finite elements, each with uniform conductivity, and Laplace's Equation is solved at all discrete elements. Creation and solving of the FEM model followed existing procedures [14,23]. Briefly, we used ScanIP to generate a mesh with tetrahedral elements using adaptive meshing to limit the size of the total model while preserving a 1 mm resolution. The FEM calculation was performed with Abaqus (Simulia, Providence, RI, USA) using automated batch-processing scripts. The standard Laplace equation was solved using conjugate gradients iterative solver. For each of the 336 HD electrodes the applied current was set to 1 mA by adjusting current-density to the exact area of each electrode. Abaqus solutions on the tetrahedral grid were read back into MATLAB and interpolated on the original 3D regular grid of the MRI segmentation.

Pad sampling and electrode weighting
For comparison, ten continuous pad configurations were also solved ( figure 1(b)). The pad-electrodes were placed manually on the scalp using ScanIP. The selection of these pad configurations, including varying pad sizes, is based on previous clinical use for the configurations C3-C4 [24], F3-F4 [4], Oz-Cz [25] as well as arbitrary combinations of other locations to quantify the possible variability in the results. Subsequently, these solutions were approximated by sampling the pad with the HD electrodes, i.e. the sampled solutions are a linear superposition of the solutions of the subset of HD electrodes overlapping in location with the padelectrode. An electrode was included in the sampling subset if at least half of its area was covered by the pad's area (figure 2). With this criterion and given the fixed locations of the sampling electrodes, a varying number of electrodes resulted for each configuration ( figure 1(b)). Both, overlap of sampling electrodes as well as their number affects the accuracy of sampling. This will be discussed in detail with figures 7 and 8. Notice also that the size of the pad electrodes for the Oz-Cz configuration differs from the other configurations. With this we can analyze the sampling technique for another pad size, i.e. a rectangular pad of 5 × 7 cm 2 , instead of only quadratic, 5 × 5 cm 2 pads. The specific choice of pad size was based on previous modeling studies that used those dimensions [25].
Simulations show that current-density distribution on a continuous pad is not uniform as shown in [26]. To explore whether emulating such an uneven distribution on the array could improve performance, we tested three different methods of distributing current among the subset of HD electrodes (uniform, non-uniform, and optimal).
Uniform. The simplest approach is an uniform weighting whereby every single electrode is set to the same current intensity value ( figure 3(a)).
Non-uniform. In a continuous pad, currents are not uniformly distributed (see for instance figure 3(e)). To emulate this behavior we tested three simplified non-even distributions of currents on the available subset of electrodes (figures 3(b)-3(d)). Electrodes were defined as corner, border or center electrodes depending on their distance from the center of the pad (corner: outside of outer ellipse, center: inside the inner ellipse, border: between the two ellipses; the axes of the outer ellipse are given by the pad dimensions and the axes of the inner ellipse are 75% from these values). The following weightings of current-density for corner-border-center were tested: 3-2-2 ( figure 3(b)), 6-3-2 (figure 3(c)) and 15-5-2 ( figure 3(d)). The specific values were selected to roughly approximate the decay we observed in current-density for different pad sizes (larger pads decay faster). The uniform distribution can be represented by the weights 1-1-1.
Optimal. Since the solutions for the continuous pad configuration are available, we can also compute a best-case solution by optimizing the weightings using least-squares optimization, which we constrain to have correct signs and to add up to a total of zero current (i.e. electrodes in a pad are either cathodal or anodal and all current entering must also exit the head). This can be implemented with linearly constrained least-squares [23]: subject to s i > 0, s j < 0 and s T · 1 = 0, (1) where E is the matrix with the electrode's solutions used to sample the pad, p is the actual pad solution arranged as a vector and s are the weighting factors for each electrode which are to be optimized. i, j indicate the index of the anodal and cathodal electrodes, respectively.

Performance evaluation
The solutions for the continuous pad and sampled pad were compared in terms of field magnitude and location of peak stimulation, which is perhaps a more relevant criterion for clinical applications. The electric fields obtained for the continuous and sampled pad are denoted by E p and E s , respectively. With these the relative error is defined as: where 2 denotes the square sum over voxels and xyz dimensions of the field vector. The calculation of relative errors was performed for gray and WM voxels separately. Areas of peak intensity (upper quartile of field intensities) were compared using the Jaccard index: P p , P s are the peak masks from pad and sampled pad results, respectively. A Jaccard index close to 1 indicates high overlap while an index close to 0 denotes no overlap.

Results
We leverage the linearity of the Laplace equation to simulate a single pad-electrode with a linear superposition of smaller electrodes. Two important questions arise in this context: how accurately do the resulting electric fields with this sampled solution replicate the fields obtained with the continuous pad-electrodes? Furthermore, how should one distribute the current among the subset of HD electrodes to achieve a faithful replication of the continuous pad-electrodes? To answer these questions we simulated ten pad configurations. The comparison of the current distribution between continuous pad and sampled pad reveals the following.

Intensity on the brain is generally accurate with an average error of 5%
The results for field magnitude on the brain are exemplified for one pad configuration in figure 4 showing a good correspondence of the sampled pad with the continuous pad. To quantify the difference we calculate the relative error (equation (2)) for all configurations, using different possible current distribution among the subset of sampling electrodes (figure 3). With a uniform current distribution the relative error in field magnitude is 5% ± 0.68% for GM and 4.6% ± 0.7% for WM (mean ± std across ten pad configurations). Other current distributions give similar results with the exception of the optimal solutions, as expected, as it uses additional information ( figure 5(a), section 3.4).
The relative error calculation for different tissue types, using the uniform distribution among the HD electrodes, reveals an increasing error for tissues that are closer to the electrodes ( figure 6(a)). This is an expected outcome: as one moves further away of the electrodes the blurring of current distributions (primarily at skin/bone and CSF/bone boundaries) reduces the effect of discrete sampling on the surface.

Accuracy in terms of the location of peak intensity on the brain is generally around 94%
To determine the locations of maximal stimulation and whether they differ significantly for the sampled solution, we measured the Jaccard index (equation (3)) comparing areas of peak field intensity (top 25%). Results indicate a very good overlap with    a Jaccard index of 94.3% ± 2% for GM and WM combined ( figure 5(b)) using the uniform distribution. The Jaccard indices for the different tissue types, also using the uniform distribution, show an increasing overlap of peak intensity locations for tissues that are closer to the electrodes ( figure 6(b)). Again, this is expected as the sampling has been selected precisely to overlap in area, while at a distance the blurring of current flow may lead to a reduced overlap.

Sampled region requires complete coverage
The ten different pad configurations are all comparable in terms of errors, which are in the range of 4%-6% (figure 7(c)). A major determinant for the adequacy of this sampling approach is a complete coverage of the pad with the sampling electrodes. An example of this problem is shown in figure 7(a) where partial coverage leads to substantially higher relative error in field magnitude (figure 7(c)) and a lower Jaccard index for peak intensity ( figure 7(d)).
To investigate the necessary resolution of HD electrodes, we computed the error rates for different number of electrodes used to sample a pad-electrode. Starting with two electrodes (one for each pad) we incrementally added electrodes and calculated the relative error and locations of peak intensity. Electrodes were selected at random and so this evaluation was repeated 100 times per configuration to obtain a mean behavior with increasing electrode counts. Evidently, increasing the number of HD electrodes leads to a decreasing relative error ( figure 8(a)) and an increasing Jaccard index, indicating a higher overlap of peak intensity areas ( figure 8(b)). One should use caution when extrapolating this data: these simulations used a constant electrode size, and in addition, as the number of electrodes increases the spatial overlap of the set of electrodes with the exact shape of the pad improves.

A non-uniform distribution of currents across the HD electrodes does not meaningfully outperform a simple uniform distribution
To determine if the different methods of distributing currents among the sample electrodes differ, we performed a repeatedmeasures one-way analysis of variances (ANOVA) with five conditions (uniform distribution, 3-2-2-, 6-3-2-, 15-5-2weighting and optimal distribution). A significance level of 5% was considered for all hypothesis tests. The test showed a significant difference (p = 2.2 × 10 −16 ) between the various methods. Obviously, the optimal solution outperforms the others, since it uses the pad solution as a prior information. Inspection of these optimal solutions reveals that they are not readily predictable and thus this performance cannot be trivially achieved in practice. Leaving out the optimal solution and conducting an ANOVA for the four remaining methods still results in a significant difference (p = 0.0275). A comparison of the mean values shows that the 15-5-2weighting produces the highest error and leads to this test result. Pairwise t-tests reveal that there is no significant difference between uniform and 3-2-2-weighting (p = 0.36), 6-3-2-weighting (p = 0.85) and 15-5-2-weighting (p = 0.17). Comparable results were obtained for the WM tissue. In summary, none of the tested weight-distributions outperform the naive uniform distribution.

Discussion
Currently the computation and preparation time needed to evaluate a specific electrode configuration is approximately 6 h for an experienced technician in a specialized laboratory and assuming segmented head models are already available. In contrast, we envision a two part process: step 1 evaluates the FEM model with high-end computers and software to precompute solutions for each HD electrode in the array. In step 2 practitioners use a thin client to quickly explore different electrode configurations.
With this approach the initial calculation of FEM solutions (step 1) would still require approximately 1 h per HD electrode. This computationally demanding step would not happen at the site of the end-user, which typically will not have access to such specialized software and hardware. Instead this computations would be performed where such expertise resides, i.e. at a university laboratory or commercial company. Additional expertise needed for this step 1 is the segmentation of the high-resolution MRI, which is not feasible in a clinical setting (radiologists do not have the time for detailed segmentation of the entire head anatomy, and automated techniques are still in development). However, step 2 can readily be implemented on a low-end computational platform at the site of the enduser. This step takes no more than 2 min on a typical desktop computer (5 s per electrode, with 20.3 ± 1.34 electrodes for the configurations tested here).
Initially we envision the use of standard head models with normal anatomy as we have done here (i.e. not individualized). At present much of the computational modeling that has been published and which is being used to guide treatment was developed on a few such examples heads [21,22]. With standard heads, the expensive step 1 is performed only once and step 2 is executed by many different clinical researchers that have varying targets and electrode configurations in mind.
Sampling provides a fast approach to simulate the current distribution in the head for the use of tES, with a relatively small decrement in accuracy. Sampling is sufficient thanks to the significant amount of blurring on the head surface and CSF. The focal spot of an HD electrode is broadly distributed at the level of cortex due to lateral current flow on skin and CSF. To be specific, for a spherical head model one can show that a point source on the skin surface results in a focal spot of approximately 27 • in width at the cortical surface [27].
Thus, we expect that the focal spots generated by two small HD electrodes that are 9 • apart (as in our current model) are already no longer resolvable at the cortical surface.
The uniform distribution over all electrodes leads to the most robust results with an error of no more than 6%. A corresponding graphical user interface that uses this sampling approach is forthcoming based on existing commercial software (HDExplore, Soterix Medical Inc., used for instance in [28]).
We have used a dense array with 300+ electrodes for obtain these results. Evidently performance would increase with increasing electrode count (figure 8) but one should use caution in extrapolating this data. Regardless, it is evident that increasing the number of electrodes provides diminishing returns and a 6% error is within the variance that is expected from different modeling choices (e.g. nonuniform conductivities [29,30], segmentation error [17], etc). So further increasing electrode counts may not be justified based on the current state of the art in current flow modeling.