Transcranial Electrical Stimulation (tES) is a non-invasive brain stimulation method used for stimulating neuronal activity, treating psychiatric disorders, and studying neuronal behavior by transmitting a constant low-intensity current pattern through a set of electrode patches attached to the scalp of the subject to modulate cortical excitability (Nitsche and Paulus, 2000). In tES, a volumetric current density in the brain is generated by injecting through the scalp a current pattern that can be described via different properties, including the number of active electrodes, their physical description (e.g., positioning, shape, permittivity, and impedance values), the applied stimulus waveform (e.g., amplitude, pulse shape, pulse width, and polarity), the number of stimulation sessions, and the time interval (Peterchev et al., 2012). Since different electrode montages result in distinct brain current flow, clinicians and researchers can adjust the montage to target or avoid specific brain regions in an application-specific manner.

An increasingly popular form of tES is the Transcranial Direct Current Stimulation (tDCS) method (Paulus, 2011; Moreno-Duarte et al., 2014; Thair et al., 2017; Reed and Cohen Kadosh, 2018). Compared to other non-invasive stimulation methods, the advantages of tDCS can be attributed to its inexpensive and approachable characteristics. Unlike the intricate machinery required for Transcranial Magnetic Stimulation (TMS) or the specialized frequency considerations in Transcranial Alternating Current Stimulation (tACS), tDCS involves a simpler setup—a direct current passed through scalp electrodes. This simplicity not only reduces the cost of equipment but also enhances portability, making tDCS more accessible for various settings, including home use. The simplicity and minimal training required contribute to its user-friendly nature enabling a broader range of individuals to utilize or participate in studies involving this method. Whereas tDCS is classically applied in a two-channel configuration (Kaufmann et al., 2021), its focality can be enhanced via multiple channels, which has motivated the introduction of advanced optimization methods for finding an optimal multi-channel montage (Fernandez-Corazza et al., 2020).

tES modeling involves constructing computational representations of the head and brain anatomy, simulating the distribution of electric fields. This process integrates factors such as electrode placement, tissue conductivity, and finite element method simulations to visualize and analyze the spatial distribution of the electric field within the brain. Generating a high-resolution forward model is critical for building an explicit patient-specific head model, determining optimal positioning of electrodes, and predicting electric field generation across the brain for specific stimulation configurations (Faria et al., 2011; Rampersad et al., 2013; Wagner et al., 2013). Using such a forward model, multi-electrode stimulation can be optimized via specifically designed mathematical methodology (Dmochowski et al., 2011; Ruffini et al., 2014; Guler et al., 2016; Wagner et al., 2016; Fernandez-Corazza et al., 2020), such as the recently developed convex optimization schemes including the Distributed Constrained Maximum Intensity (D-CMI) (Khan et al., 2022), and the metaheuristic L1-norm regularized L1-norm fitting (L1L1) (Galaz Prieto et al., 2022) which aim at an individualized distributional fit for a given target activity.

In this study, we aim to broaden the applicability of the linear programming (LP)-based L1L1 method for finding tES electrode montages computationally in a comprehensive manner, i.e., by evaluating the metaheuristic results and total computing time through different mathematical optimization algorithms and packages; this includes Interior-Point (IP) (Mehrotra, 1992), Primal-Simplex (PS), and Dual-Simplex (DS) (Boyd and Vandenberghe, 2004) as alternative LP algorithms, and Matlab (R2020b) from MathWorks (Zhang, 1999), MOSEK Optimization Suite (Release 9.0) (Mosek, 2019), Gurobi Optimization (9.5.1) (Gurobi Optimization LLC, 2022), SDPT3 (4.0) (Tütüncü et al., 2003), and SeDuMi (1.3.5) (Sturm, 1999; Frenk et al., 2000; Polik et al., 2007) as alternative packages. The latter two open-source alternatives are available in the CVX optimization toolbox (Grant and Boyd, 2014). We also investigate the metaheuristic hyperparameter optimization (HPO) task of L1L1 via exhaustive search (Bianchi et al., 2009) and recursive search (Je and Park, 2013) with heuristic direct search as a reference algorithm (Bogani et al., 2009).

Our results suggest that the performance differences between the above-mentioned optimization packages, algorithms, and search methodology can be crucial regarding the optimization results, focality stimulation current, and the availability of active channels in the montage. Moreover, exhaustive and recursive search methods can also be considered preferable to heuristic direct search in terms of their overall reliability and predictability.

In tES, a real L × 1 current pattern y is injected into the subject's head through a set of contact electrodes attached to the scalp. These electrodes, ranging from 0.5 to 4.0 milliamperes (mA) (Zaghi et al., 2010; Khadka et al., 2020; Workman et al., 2020), form what is known as an electrode montage and are responsible for distributing the injected volumetric current density–measured in ampere per square meter (A/m²)–throughout the scalp, skull, cerebrospinal fluid (CSF), and brain components, including cortical and subcortical brain structures. The governing linear system is of the form

where L^ is a real N × L lead field matrix (forward mapping) that describes the relationship between the y, and x^ is a real N × 1 discretized volume current density vector. The linear system (Equation 1) is re-interpreted component-wise as the focused field L^1y=x^1, where the target field has non-zero values, and the nuisance field L^2y=0, where it vanishes. Detailed mathematical definition of the lead field matrix refer to Appendix A. Forward model in Galaz Prieto et al. (2022).

The optimization problem needs to find the best matching between y, and the focused field via Ly = x, where the projection of the focused field into the direction of the target constitutes the first component as

with P denoting a matrix that projects a vector into the direction of x^1. The target amplitude ||x₁||₂ is set as 3.85 A/m² which is an approximation of the excitation current threshold for nerve fibers of the upper limb area of the motor cortex (Kowalski et al., 2002).

The exhaustive search proved to be a reliable method for experimental benchmarking when the required tES montage requires careful design for clinical applications. By presenting the exhaustive search results in the form of a heatmap with a coarse grid of κ = 15 (Figure 3A), we can pinpoint the (α, ε) region where the focused current amplitude reaches its maximum. Despite a significantly increased number of evaluations, with a finer grid of κ = 40 (Figure 3B), we can further determine a more detailed optimal area. This area corresponds to the Cartesian product of α_m ranging from −71 to −50 dB and ε_n from 0 to −98 dB. In this context, a high current injection montage, denoted by Γ_max (yellow star), is positioned at the peak of the amplitude, while focality-based montages, Θ_max (purple star), adhere closely. However, these focality-based montages are slightly deviated due to the influence of the nuisance field, despite being relatively close, as determined by a threshold condition corresponding to 75% of the maximum amplitude achievable with the two-patch bipolar tES montage. In comparison between these grid resolutions, one can observe slight enhancements in amplitude, increased optimization accuracy, and improved numerical stability in the latter case. These aspects are far more noticeable with Dual- and Primal-Simplex methods than with the Interior-Point, which yields overall smoother results with fewer drastic deviations.

Figure 4 delineates the performance nuances among optimization strategies. The whiskers along the stems signify a second-order Taylor's polynomial estimate, reflecting the maximum deviation within half lattice units distance from the optimizer. The reciprocity principle reference for Γ_max is represented by a horizontal black dashed line, and the number of non-zero (NNZ) channels required for a tES montage is depicted on the right side of each corresponding stem. The solvers are sorted in ascending order based on their performance, with the exhaustive search κ = 40 grid (blue) serving as the point of reference. Both the direct and recursive search techniques adeptly uncover optimal (α, ε) solutions for Θ_max and Γ_max, yielding a substantial reduction in total computing time compared to the specified hyperparameter space.

Due to the heuristic nature of the direct search, and to assert the efficacy of the said technique, we performed a series of trial runs by setting the initial point to the center of the search space. In these trials, the number of objective function evaluations varied, ranging from 25 to 54 trials, with a 33.8 mean among the evaluations (see Table 2). While the number of function evaluations was slightly higher, the quality of the results was nearly on par with those obtained with the recursive search with a search window of K~=3. With IP solvers, the search runs were mostly successful, while PS and DS tended to fail to find a feasible optimizer candidate.

Due to its relatively fast performance among interior-point methods, we applied MOSEK IP to evaluate topographical maps of stimulus focality Θ_max (Figure 5A) and current density Γ_max (Figure 5B) for an exhaustive search κ = 15 and a recursive search K~=3. Overall, the results of the recursion were close to the outcome of the exhaustive search. Thus, the topographical differences between the different approaches of this study were observed to be minor.

By limiting the search space only to a narrower subspace is a simple countermeasure for dealing with the disadvantages of the exhaustive search. With a κ = 15 grid as a reference, it can take approximately 850 seconds to perform a complete search for the first stage, while the second stage only takes roughly 15% of that time since it uses a limited lead field following from the limited number of active electrodes. As an alternative approach, the direct and recursive search seemed to perform well compared to the number of objective function evaluations made during the search process (Figure 6). In particular, MOSEK turned out to be the superior choice, with MOSEK DS being the fastest one. The computing time for Gurobi IP was close to that of MOSEK IP, and Gurobi DS, PS, and Matlab IPL and DS required approximately three times the time. The slowest-performing SDPT3 and SeDuMi took as much as six times the run time of MOSEK IP. Overall, the simplex methods applied to the L1L1 optimization scheme deliver faster yet less accurate solutions than Interior-Point (IP) for focality-based montages, while minor differences can be found for intensity-based solutions.

In this study, we analyzed the numerical and computational performance of exhaustive search, direct search, and recursive search techniques to find an optimal stimulation focality Θ and current density Γ for solving the L1L1 optimization problem for non-invasive transcranial electrical stimulation (tES) current injection. This analysis was motivated by our earlier results in (Galaz Prieto et al., 2022) which suggested that the L1L1 method provides a theoretically attractive approach for obtaining a high-gain focal stimulus as compared to complex L2-norm fitting and regularized least squares techniques (Dmochowski et al., 2011; Wagner et al., 2016).

The reciprocity principle, as outlined by Fernandez-Corazza et al. (2020), served as a reference technique. Its validity was shown for the present tES lead field matrix L (see Section 2.3.1). When focusing on a specific target region, the current injection pattern from a two-patch tES montage aligns with the maximum intensity achievable through this principle. Essentially, this involves selecting the two electrodes with the highest absolute back-projected currents. With the absence of nuisance field constraints, the L1L1 solution was observed to agree with the reciprocity principle if the aforementioned algorithmic aspects were handled appropriately.

Decisive aspects for a successful outcome of L1L1 were found to be the choice of the optimization package, algorithm, and search routine, which significantly affect both the performance of the metaheuristic optimization process and output. To enlighten this aspect, we covered the performance of several Interior-Point (IP) (Mehrotra, 1992), Dual-Simplex (DS), and Primal-Simplex (PS) (Boyd and Vandenberghe, 2004) methods from different open-source and commercial optimization toolboxes. We tested the L1L1 method using the commercial solvers of MOSEK Optimization Suite (Release 9) (Mosek, 2019) and Gurobi Optimization (9.5.1) (Gurobi Optimization LLC, 2022), and compared them to the open-source alternatives (Grant and Boyd, 2014) SDPT3 (4.0) (Tütüncü et al., 2003) and SeDuMi (1.3.5) (Sturm, 1999; Frenk et al., 2000; Polik et al., 2007) as well as Matlab R2020b's (MathWorks) Interior-Point-Legacy (IPL) algorithm, which originates from the open LIPSOL (Zhang, 1999) toolbox. We selected the IPL algorithm since we experienced stagnation with Matlab's main IP algorithm, which did not return any appropriate results.

Based on the results, we consider Gurobi IP to be the preferable choice in both optimization stages, considering Θ_max and Γ_max in each tested target region of interest and, as it was also overall the fastest of the IP solvers. While the best-performing solvers show that the L1L1 method is suitable for maximizing focality and intensity, a few did not find the bipolar current pattern that maximizes Γ_max. Notably, SDTP3 did not find a bipolar pattern at all, verifying our earlier hypothesis (Galaz Prieto et al., 2022) that the performance of L1L1 might be highly solver-based. Part of the discrepancies between the optimization methods can be explained by a different sensitivity with respect to parameter variation or the resolution of the lattice.

From a computational complexity standpoint, the exhaustive search method can be applied for benchmarking purposes. In contrast, a recursive search proves an advantageous alternative and is competitively on par with the direct search technique, each one applied in this study. This equivalence arises from both methods converging toward the most suitable regularization parameter α and nuisance threshold ε values in a comparably controlled manner. Notably, the computational complexity of recursive search remains consistent across various optimization runs, in contrast to the variability observed in the direct search. Results comparable to those obtained through exhaustive search can be attained with a reduced-resolution search window of, say, size K~=3, representing a substantial acceleration in comparison to exhaustive search. Furthermore, the recursive approach demonstrates both numerical stability and convergence towards exhaustive search results, both at individual data points and in the overall topographical context, as the probing lattice size increases.

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

FG: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing—original draft. MS: Formal analysis, Project administration, Supervision, Validation, Software, Visualization, Writing—original draft, Writing—review & editing. SP: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing—review & editing.