# An optimization approach for well-targeted transcranial direct current stimulation

###### Abstract

Transcranial direct current stimulation is a non-invasive brain stimulation technique which modifies neural excitability by providing weak currents through scalp electrodes. The aim of this study is to introduce and analyze a novel optimization method for safe and well-targeted multi-array tDCS. For optimization, we consider an optimal control problem for a Laplace equation with Neumann boundary conditions with control and point-wise gradient state constraints. We prove well-posedness results for the proposed methods and provide computer simulation results in a highly realistic six-compartment geometry-adapted hexahedral head model. For discretization of the proposed minimization problem the finite element method is employed and the existence of at least one minimizer to the discretized optimization problem is shown. For numerical solution of the corresponding discretized problem we employ the alternating direction method of multipliers and comprehensively examine the cortical current flow field with regard to focality, target intensity and orientation. The numerical results reveal that the optimized current flow fields show significantly higher focality and, in most cases, higher directional agreement to the target vector in comparison to standard bipolar electrode montages.

Key words. transcranial direct current stimulation, optimization, Neumann boundary control, gradient state constraint, sparsity optimal control, existence

AMS subject classifications.

## 1 Introduction

Transcranial direct current stimulation (tDCS) is a non-invasive, inexpensive and easy-to-perform brain stimulation technique which modifies neural excitability [23].
Changes in neural membrane potentials are induced in a polarity-dependent manner, for example, in a motor cortex study, anodal stimulation right over the motor cortex enhances cortical excitability, whereas cathodal stimulation inhibits it [23]. Recently, tDCS has been applied successfully in the treatment of neurological and neuropsychiatric disorders such as epilepsy [14], depression [6] and Alzheimer disease [13]. The effects of tDCS can be preserved for more than one hour after stimulation [24].

The conventional strategy is to apply the current density via two large electrodes, with the active electrode (anode) to be placed above the presumed
target region and the reference electrode (cathode) far away from the target region [23, 24]. Accurate and detailed finite element (FE) head models have been created to investigate the induced current density distribution [32, 10, 30]. While significant effects of stimulation as compared to sham were reported [6, 14, 23], computer simulation studies have revealed that the induced cortical current flow fields are very widespread with often strongest current density amplitudes in non-target brain regions [32, 30]. It is therefore a matter of debate whether the effects of stimulation are driven by the target brain region or elicited by adjacent cortical lobes [30].

In order to overcome the limitations of conventional bipolar electrode montages, algorithmic-based sensor optimization approaches were presented [10, 30, 21, 28]. Im and colleagues [21] searched for two electrode locations which generate maximal current flow towards a certain target direction. Ruffini and colleagues [28] described a method for optimizing the configuration of multifocal tCS for stimulation of brain networks, represented by spatially extended cortical targets.
Dmochowski and colleagues [10] used a multi-channel array consisting of 64 fixed electrode locations to calculate optimized stimulation protocols for presumed target regions. They employed radial and tangential targets and reported that compared with conventional electrode montages their optimization approach achieved electric fields which exhibit simultaneously greater focality and target intensity at cortical targets using the same total current applied.

While the optimization results are very promising and a pilot study using multi-array tDCS devices for rehabilitation after stroke has been conducted [11], to our knowledge, there is currently no study providing an in-depth analysis of mathematical models for multi-array tDCS optimization.

Given a volume conductor model with a fixed electrode arrangement and a target vector in with being the target cortical area, the optimization approach estimates an optimal applied current pattern at the fixed electrodes. A Laplace equation with inhomogeneous Neumann boundary conditions is used to calculate the induced current density distribution [32] which is to be controlled by the boundary condition to ensure safe and focused stimulation. Therefore, the optimization problem for tDCS is in the class of control problems with Neumann boundary conditions [22].

In this paper, we modeled the electrodes with the point electrode model (PEM)
[26, 1, 2]. They can, however, also be modeled using a complete electrode model (CEM) [26, 9, 12, 1, 2]. In [1, 2], it has been shown, that CEM and PEM only lead to small differences which are mainly situated locally around the electrodes and are very small in the brain region. Based on these results, the application of PEM
is expected to result in negligible differences to the CEM and should thus provide a sufficiently accurate modeling of the current density within the brain region.

This paper is organized as follows: In Section 2 we establish the mathematical model of tDCS and introduce the optimization problem for multi-channel tDCS. In Section 3 we show existence of at least one minimizer to the optimization problem and the FE method is used for numerical discretization. In Section 4 we derive the algorithm for sensor optimization. Section 5 provides optimization results in a highly-realistic six-compartment head model (skin, skull compacta, skull spongiosa, CSF, gray and white matter) with white matter anisotropy. Furthermore, a conclusion and outlook section is presented.

## 2 Mathematical Model of tDCS

In order to calculate the current flow field induced by tDCS, the quasistatic approximation to Maxwell’s equations is applied [32]. This yields the tDCS forward problem

with being the electric potential, being an anisotropic conductivity tensor, being the applied current pattern at the electrodes with non-zero values only at the electrode surfaces, being the outward normal vector and being a part of the boundary of the domain . Furthermore, a Dirichlet boundary condition on the remaining part is used to ensure that a solution to the tDCS forward problem is unique. We assume all parts of the boundary to have nonzero measure and to be of Lipschitz regularity. We introduce the Sobolev space

with being the standard Sobolev space for Neumann boundary values on . The precise definition of is to be the dual space of consisting of the Dirichlet traces of -functions in . The integral in the above definition has to be interpreted as a duality product with the constant function equal to one, which is in . The weak formulation of the forward problem is given by

for all test functions with vanishing Dirichlet trace on . Note that the boundary integral on the right-hand side is only defined if , in general it is to be replaced by the duality product between and

Under standard and naturally satisfied regularity assumptions, there exists a unique solution to the tDCS forward problem, which can be shown by the Lax-Milgram Lemma:

###### Theorem 2.1

Let with and let . Then there exists a solution to the tDCS forward problem. The solution is unique if has positive measure. Moreover, the following norm estimate holds:

For current density optimization, a multi-channel array consisting of 74 fixed electrode locations (the locations of an extended 10-10 EEG electrode system) is used and optimized applied current patterns are calculated for presumed targets. Figure LABEL:Fig0 illustrates an overview of the optimization setup in a two-dimensional model. Optimally, the induced current density distribution is maximal in the target region and zero in non-target regions . Physically, a focal stimulation without stimulating non-target regions in not possible as the current has to flow through the volume to reach the target region. Therefore, the current density is restricted by such that

Secondly, the total current applied to all electrodes is limited to 2 mA, a commonly used safety criterion [10].

For effective stimulation, the optimized current density distribution should be oriented perpendicularly to the cortical surface [5, 8]. Besides the correct target location, also the target direction is important as shown by Bindman and colleagues [5] and Creutzfeldt and colleagues [8] who were able to demonstrate in physical measurements in rats and cats, respectively, that the neural firing rate is strongly influenced by the direction of the field to the cortical surface. Perpendicularly-inwards (and parallel to the long apical dendrites of the large pyramidal cells in cortical layer V) stimulation, i.e. anodal stimulation, strongly enhanced the activity of the cortical neurons, whereas perpendicularly-outwards stimulation, i.e. cathodal stimulation, inhibited it. We thus maximize with being a unit vector perpendicularly-inwards-oriented to the cortical surface [5, 8].

In the following we reformulate the current density constraint in the whole domain as
with a weight in and in . Thus, we consider the
constrained optimization problem

subject to | ||||

which is a control problem with Neumann boundary conditions [22]. Note that safety limitation on the inflow current is a control constraint, while the limit on the current density outside the target region is a state constraint. Also, the state constraint with ensures that is bounded in and due to the embedding theorem also bounded in . This implies that and due to the trace theorem . Since this implies that is bounded in .

###### Remark 2.2

The constrained optimization problem is a rather challenging optimal control problem. Firstly, the current density distribution is optimized, while many other applications require to optimize the potential field . Secondly, a point-wise gradient state constraint is used. While optimal control problems with state constraints were frequently used [17], not much attention was given to gradient state constraints only until recently [29, 25, 33] focusing on distributed control in the domain however. Thirdly, we look for a control function vanishing on large parts of the boundary, known as sparsity optimal control problem [19]. Indeed, we rather look for as a combination of concentrated measures rather than an -function, so the above integral formulation has to be understood as formal notation for the total variation of the Radon measure identified with . This will be made precise in the next section.

Applying Gauss’ Theorem to provides a further motivation for the proposed objective functional. Assuming and locally constant we have

In order to maximize the current densities along the target direction, the potential should be maximized where the target vector is oriented parallel to the outward normal vector on and minimized where it is antiparallel.

## 3 Optimal Control Problem

The goal of this section is to provide insight into the existence theory of at least one minimizer to the tDCS optimization problem . Furthermore, we show existence of at least one minimizer to a simplified optimization problem without control constraint at the electrodes. Finally, the finite element method is used for numerical discretization of the considered optimization procedure.

### 3.1 Continuous formulation

We now reconsider the optimization problem . In order to provide a rigorous formulation of the safety constraint we interpret as a Radon measure in the space , whose norm, i.e. the total variation of a measure, is an appropriate replacement for the -norm. The constraint then becomes

In order to obtain a unified formulation also allowing for other design constraints respectively a discrete electrode setup, we introduce a feasible set . To further simplify notation we define an operator via

Since , the boundary value problem for the Poisson equation has a unique solution bounded by a multiple of the norm of as guaranteed by Theorem LABEL:thmwelldef. This implies that is a well-defined linear operator. Hence we can substitute in , which leads to the following equivalent reformulation

subject to | ||||

The existence of at least one minimizer to is the topic of the following theorem.

###### Lemma 3.1

Let be bounded away from zero in and let be such that satisfies the constraint almost everywhere in . Then there exists a constant depending only on the given data , , , and , such that

\hb@xt@.01(3.1) |

Proof. We have

By the trace theorem we can write each such as the Dirichlet trace of some with vanishing trace on and for a universal constant in the trace theorem depending only on and . Hence, by the weak formulation of the elliptic partial differential equation satisfied by with , we find

Using the Cauchy-Schwarz inequality and an elementary estimate of the -norm of in terms of its supremum norm, we find

and we see that only depends on , , , and .

###### Theorem 3.2

Let be weakly closed in , let and let . Then there exists at least one minimizer of the convex optimization problem .

Proof. An easy computation reveals that is feasible for , i.e. the set of feasible points is not empty. Moreover, it is closed in the weak-star topology of , which can be seen as follows: First of all the norm in is lower semicontinuous with respect to weak-star convergence in this space, which implies that a limit of feasible points satisfies the safety constraint as well. Moreover, weak convergence of a sequence in implies weak convergence of in , which implies that the limit satisfies the state constraint due to the weak closedness of pointwise constraints in . Moreover, the feasible set is obviously bounded in due to the safety constraint and in due to Lemma LABEL:newlemma. Hence, the Banach-Alaoglu Theorem implies compactness in the weak-star topology. Since the objective is a bounded linear functional and in particular weak-star lower semicontinuous we obtain the existence of a minimizer by standard reasoning.

We mention that our existence proof implicitely uses additional regularity induced by the gradient state constraint as done also in [33], however we do not discuss further regularity issues of the solution, which is more involved and possibly not to be expected with the nonsmooth terms in the objective functional.

We observe that the objective as well as all terms involved in the constraints of the optimal control problem are one-homogeneous. This allows us to renormalize the unknown and in this way eliminate the safety constraint. This motivates the following simplified version:

subject to |

By renormalizing a minimizer of we can obtain a minimizer for the full optimization problem with safety constraint for the control, which of course only holds if the total variation of the minimizer is finite:

###### Theorem 3.3

Let be a threshold parameter and be a minimizer of the simplified minimization problem . Then is a minimizer of with .

Proof. An easy calculation reveals that fulfills the control and state constraint of with . Assume that is not a minimizer of . Then there exists with smaller functional value, and
is feasible for with lower functional value than , which contradicts the optimality of the latter.

At a first glance the usefulness of Theorem LABEL:PeQe might appear limited, since the relation between the constraint bounds and is quite implicit, depending on the minimizer of the second problem. However, since there is no natural choice of the bound in either case one can choose or in appropriate ranges equally well. It turns out however that is somewhat easier in particular with respect to numerical solutions. Moreover, it offers a better position to introduce sparsity constraints on the control by additional penalization, which is the subject of the following discussion.

### 3.2 Penalized Problem Formulation

In order to control the applied current pattern at the fixed electrodes, the minimization problem is extended by two penalties. While an term is introduced to penalize the energy of the applied current pattern, an penalty is used to minimize the number of active electrodes in the minimization procedure. The latter has to be reinterpreted in the same way as the safety constraint in terms of Radon measures. This leads to the following minimization problem

subject to |

We call the minimization problem with and to be the regularized optimization procedure (L2R) and the regularized optimization procedure (L1R), respectively. If both penalties are active, the penalty term is also known as elastic net.

Since the penalty adds strict convexity to the problem, we can obtain uniqueness of a minimizer together with existence via analogous arguments as above:

###### Theorem 3.4

Let be weakly closed in , , and let and . Then there exists a unique minimizer to the L2R constrained optimization problem . In the latter statement, can even be generalized to .

For the problem with pure -penalization the proof of existence is almost identical to the one of Theorem LABEL:PeQe and we hence conclude:

###### Theorem 3.5

Let be weakly closed in and let and . Then there exists at least one minimizer to the L1R constrained optimization problem .

### 3.3 Discrete Problem

For numerical discretization of the considered optimization problem we use the finite element method. The solution of the forward problem is approximated by looking for with the being edge-based finite element basis functions vanishing on .

for all in the span of . Using standard techniques this allows to verify the existence and uniqueness of the discrete solution for arbitrary and to introduce a discretized solution operator

Let , denote the volume elements in the finite element discretization (tetrahedra or cuboids) and denote for a function the local mean value by

Finally we discretize with being local basis functions on . From the fact that has mean zero we can eliminate . Now we introduce mappings and incorporating the discretizations:

and define the discretized operator , represented by a matrix in .

As the target vector is only defined in the target region, we define a continuous mapping :

and replace by , since = 0 in . We now introduce the discretized optimization problem as

subject to |

where is a constant approximation of in (in our examples defined piecewise constant anyway).

The existence (and potential uniqueness) of minimizers can now be verified exactly as for the continuous case above, defining an appropriate discrete version of the -norm on

e.g. via

For the sake of brevity we do not discuss the issue of, which can, e.g., be shown via -convergence of the functionals involved.

## 4 Numerical Optimization

For numerical solution of the corresponding discretized problem we employ the alternating direction method of multipliers (ADMM) [7]. The ADMM is a variant of the
Augmented Lagrangian method, a class for solving constrained optimization problems. The method combines important convergence properties (no strict convergence or finiteness is required) and the decomposability of the dual ascent method [7].

To solve the discretized optimization problem , we substitute and and obtain the following Lagrangian which is to be minimized

with and being the augmented Lagrangian parameters and the dual variables, respectively [7].

-step: Rearranging and neglecting terms without , putting the integrals together and expanding leads to a quadratic minimization problem

Differentiating with respect to and setting the equation system to be zero results in

-step: Rearranging and neglecting terms and putting the integrals together results in

which can be solved analytically as follows

-step: Rearranging the terms, leaving out the terms without and putting the integrals together leads to

The solution to this equation is given as

Dual update: Finally, the dual variables are updated

We mention that a complete convergence analysis can be carried out following the arguments in [7]. Algorithm LABEL:AlgotCSOpt depicts an overview of the optimization steps for the numerical solution of the discretized problem . Note that the Euclidean norm was used for the stopping criterion in Step 3.

## 5 Results and Discussion

0.0001 | 0.0005 | 0.001 | 0.005 | 0.01 | |
---|---|---|---|---|---|

0.00332 | 0.0085 | 0.0154 | 0.0175 | 0.0208 | |

Iterations | 132 | 53 | 33 | 203 | 143 |

120.15 | 233.96 | 352.279 | 1312.2 | 1919.93 |

We generated a highly realistic geometry-adapted six-compartment (skin, skull compacta, skull spongiosa, CSF, gray and white matter) finite element head model from a T1- and a
T2-weighted magnetic resonance image (MRI). For the compartments skin, skull compacta, skull spongiosa, CSF and gray matter we employed conductivity values = , , , and , respectively [18, 27, 32]. For the modeling of the white matter anisotropy, a diffusion tensor MRI was used and the effective medium approach was applied. The effective medium approach assumes a linear relationship between the effective electrical conductivity tensor and the effective water diffusion tensor in white matter [31], resulting in a mean conductivity value for the white matter compartment [32]. As the current density amplitudes in the skin and CSF compartments are much stronger when compared to the brain compartments [32], we only visualized the current densities in the brain to enable best visibility of the cortical current flow pattern. The software package SCIRun was used for visualization [9].

Using this highly realistic volume conductor model we comprehensively examine the optimized current densities in the brain with regard to focality, target intensity and orientation. Firstly, the
maximal current density in non-target regions is investigated for the discretized optimization problem and four target vectors are used to calculate optimized stimulation protocols. For all simulations, the target area always consists of the elements corresponding to the target vectors and only the volume conductor elements in the brain are used for optimization, i.e., only in the brain compartments and in the target region and in the CSF, skin and skull compartments.

In order to investigate the current flow field that is induced by an optimized standard bipolar electrode montage, we use a maximum two electrode (M2E) approach. This approach stimulates only the main positive (anode) and the main negative electrode (cathode) of the L1R optimized stimulation protocols with a total current of . In order to quantify the optimized current flow fields, we calculate the current densities in the direction of the target vectors () and the percentage of current density that is oriented parallel to the target vector (), as shown in columns 4 and 5 in Table LABEL:Quantification, respectively.

### 5.1 Maximal current density in non-target regions

In this section, the L1R approximated discretized optimization problem with = 0.001 is used to calculate optimized stimulation protocols for a set of 924 tangential (parallel to the inner skull surface) and radial (perpendicular to the inner skull surface) target vectors and Theorem LABEL:PeQe is applied to estimate the averaged maximal current density in non-target regions as shown in Tables LABEL:TabEpsTang and LABEL:TabEpsRad. Furthermore, the number of iterations until a minimum was found and the scaling factor is depicted. As can be seen in the tables, for the threshold value = 0.001, the safety value ensures that the maximal current density in the brain compartment is not dangerous for the subject. Furthermore, the number of iterations until a minimum was found is lowest for = 0.001. For this reason, a threshold value of = 0.001 is thus used in this study as it allows safe and well-targeted stimulation and fast and robust computation of the stimulation protocol. Note that the number of iterations until a minimum is found increases with decreasing mesh size. However, because the maximal resolution of state-of-the-art whole-head MRI sequences is 1mm, we did not investigate mesh sizes smaller than 1mm.

0.0001 | 0.0005 | 0.001 | 0.005 | 0.01 | |
---|---|---|---|---|---|

0.0043 | 0.0114 | 0.015 | 0.023 | 0.025 | |

Iterations | 134 | 36 | 35 | 68 | 118 |

92.03 | 176.08 | 267.06 | 887.08 | 1597.73 |

### 5.2 Mainly tangential target vector

Figure LABEL:FigTangOpt displays the optimized current densities (for better visibility, we used a different scaling for the different rows) and the corresponding stimulation protocols (same scaling) for a mainly tangential target vector as shown in Figures LABEL:FigTangOpt-A1 and -A2. As can be seen in Figures LABEL:FigTangOpt-B1 and LABEL:FigTangOpt-C1, the optimized current flow fields using the L2R and L1R approaches, respectively, show high focality and maximal current densities can be observed in the target region. While the L2R and L1R approaches lead to target current densities of and (Table LABEL:Quantification, second column), current densities in non-target brain regions are rather weak (Figs.LABEL:FigTangOpt-B1 and -C1). Overall, due to the rather widespread applied current pattern at the fixed electrodes (Fig.LABEL:FigTangOpt-B2), the L2R optimized current flow field (Fig.LABEL:FigTangOpt-B1) has smaller amplitude when compared to the L1R one (Fig.LABEL:FigTangOpt-C1). On the other hand, the L1R stimulation protocol (Fig.LABEL:FigTangOpt-C2) shows high focality with mainly two active electrodes, while only very weak compensating currents are injected at the neighboring electrodes. The M2E approach results in higher target current densities of (Table LABEL:Quantification, second column). However, using the M2E approach, relatively strong current densities can also be noted in non-target regions, especially in pyramidal tract and deeper white matter regions (Fig.LABEL:FigTangOpt-D1).

With a value of 86.4, 86.8 and 85.9 for L2R, L1R and M2E (Table LABEL:Quantification, fifth column), the current densities of all three approaches are mainly oriented parallel to the target vector, leading to values, i.e., current densities along the target direction, of nd , resp. (Table LABEL:Quantification, fourth column). While the values between L1R and M2E are thus less than a factor of 2 different, the averaged current density amplitudes in non-target regions is about 5.3 times higher when using the M2E approach (Table LABEL:Quantification, third column). The L1R optimized current flow field thus shows significantly higher focality and also slightly better parallelity to the target vector in comparison to the bipolar electrode montage M2E. However, if no multi-channel tDCS stimulator is available, the M2E approach provides an optimized bipolar electrode montage for a mainly tangential target vector.

Target | ||||
---|---|---|---|---|

tangential L2R | 0.022 | 0.00144 | 0.019 | 86.4 |

tangential L1R | 0.038 | 0.00151 | 0.033 | 86.8 |

tangential M2E | 0.071 | 0.0080 | 0.061 | 85.9 |

radial L2R | 0.026 | 0.00064 | 0.025 | 96.1 |

radial L1R | 0.045 | 0.00071 | 0.043 | 95.5 |

radial M2E | 0.063 | 0.0074 | 0.048 | 76.2 |

patch L2R | 0.025 | 0.00147 | 0.022 | 88.0 |

patch L1R | 0.037 | 0.00157 | 0.033 | 89.1 |

patch M2E | 0.071 | 0.0080 | 0.062 | 87.3 |

deep L2R | 0.015 | 0.00345 | 0.013 | 86.7 |

deep L1R | 0.019 | 0.00249 | 0.018 | 94.7 |

deep M2E | 0.052 | 0.01477 | 0.049 | 94.2 |

### 5.3 Mainly radial target vector

The mainly radially oriented target is shown in Figures LABEL:FigRadOpt-A1 and -A2. Figures LABEL:FigRadOpt-B1 and LABEL:FigRadOpt-C1 depict the optimized current densities for the L2R and L1R approaches, resp.. As shown in Table LABEL:Quantification (Column 2), with a value of , the target current density for the L1R approach is more than a factor of 1.7 times stronger than with the L2R approach (). Non-target regions show only weak current densities (Figures LABEL:FigRadOpt-B1 and -C1). With a value of , the M2E approach yields the largest target intensity (Table LABEL:Quantification, second column). However, for this approach, the maximal current density in the brain does not occur in the target region, but on a neighboring gyrus in between the stimulating electrodes and in mainly tangential direction (Figure LABEL:FigRadOpt-D1). Moreover, the M2E optimized current density shows again an overall much lower focality than the L2R and L1R current flow fields.

The L2R stimulation protocol (Figure LABEL:FigRadOpt-B2) consists of an anode above the target, surrounded by four cathodes and a ring of very weak positive currents at the second neighboring electrodes, a distribution which might be best described by a sinc-function. As can be seen in Figure LABEL:FigRadOpt-C2, the L1R stimulation protocol is mainly composed of a main positive current at the electrode above the target region and four return currents applied to the surrounding electrodes.

The L2R and L1R optimized current flow fields show high directional agreement with the target vector ( value above ), with the L2R slightly outperforming the L1R approach (Table LABEL:Quantification, fifth column). With a value of only the current flow field in the M2E model shows much less directional agreement with the target vector. This results in values of 0.025, 0.043 and for the L2R, L1R and M2E approaches, resp.. In order to obtain higher values for the M2E approach, i.e., higher current densities along the target direction, the distance between anode and cathode might be enlarged, in line with Dmochowski and colleagues [10] who reported that the optimal bipolar electrode configuration for a radial target vector consists of an electrode placed directly over the target with a distant return electrode. Another interesting bipolar electrode arrangement for a radial target vector might consist of an anode over the target and a cathodal ring around the anode.

### 5.4 Extended mainly tangential target vector region

In this section, an extended target area of is used for current density optimization. The target area is centered around the location of the superficial target vector that was also used in Section LABEL:tangsingle and the corresponding target vectors are selected to be mainly tangentially oriented (Figures LABEL:FigPatch-A1 and LABEL:FigPatch-A2). As can be seen in Figures LABEL:FigPatch-B1 and LABEL:FigPatch-C1, the optimized current flow fields of the L2R and L1R approaches show high focality. L2R and L1R yield an averaged target intensity of and , resp. (Table LABEL:Quantification, second column). Current density is not restricted to but focused to the target area, especially for the L2R approach (Fig. LABEL:FigPatch-B1). In comparison to the superficial tangential target vector of Section LABEL:tangsingle, the L1R optimized averaged current density in the target area is decreased by about (from to ) and the optimized stimulation protocols are also very similar. Because the main two electrodes are taken from the L1R optimization, this implies that also the M2E stimulation protocol remains constant to the stimulation protocol from Section LABEL:tangsingle. In this way, the M2E approach yields an averaged target intensity of (Table LABEL:Quantification, second column).

With averaged values of 88.0, 89.1 and 87.3 (Table LABEL:Quantification, fifth column), the current densities are mainly oriented parallel to the target vectors with L1R performing best. The averaged current flow field intensity along the target direction is 0.022, 0.033 and for the L2R, L1R and M2E approaches, resp. (Table LABEL:Quantification, fourth column). However, for M2E also non-target regions reach significant current densities, as clearly shown in Figure LABEL:FigPatch-D1. The averaged current density amplitudes in non-target regions is 0.00157 and for the L1R and M2E approaches (Table LABEL:Quantification, third column). The L1R optimized current flow field thus shows a factor of 5.1 higher focality in comparison to bipolar electrode montage M2E. However, the M2E approach provides an optimized bipolar electrode montage for an extended target area of tangential target vectors.

### 5.5 Deep and tangential target vector

In the last simulation scenario, we investigate optimization for a deeper and mainly tangentially oriented target vector as shown in Figures LABEL:FigDeep-A1 and -A2. Figures LABEL:FigDeep-B1 and -C1 depict the optimized current density distributions when using L2R and L1R for optimization, resp.. For those approaches, the target current densities are 0.015 and , resp.. With a value of , which is more than 2.7 times the L1R value, the largest target intensity is, however, achieved with the M2E approach (Table LABEL:Quantification, second column).
values of 0.013, 0.018 and lead to values of 86.7, 94.7 and 94.2 for the L2R, L1R and M2E approaches (Table LABEL:Quantification, fourth and fifth column). The target current densities are thus for all three approaches oriented mainly parallel to the target vector.

The L2R and L1R stimulation protocols show high focality with mainly two active electrodes, while only weak compensating currents are injected at the neighboring electrodes (Figures LABEL:FigTangOpt-B2 and -C2). Similar to Section LABEL:tangsingle, the compensating currents are stronger when using the L2R optimization procedure, leading to weaker target brain current densities as compared to the L1R stimulation. In order to enable current density to penetrate into deeper brain regions, the distance between the two main stimulating electrodes is larger when compared to the superficial mainly tangential target vector from Section LABEL:tangsingle, i.e., the electrode above the target region is not used for stimulation, while a more distant electrode is used as anode.

For all three approaches, strongest current density amplitudes in the brain compartment always occur at the CSF/brain boundary above the target region (Figures LABEL:FigDeep-B1,-C1 and -D1). This is due to the fact that the potential field satisfies the maximum principle for harmonic functions which states that a non-constant function always attains its maximum at the boundary of the domain [15, Theorem 14.1]. Nevertheless, in average over all non-target regions, with a value of for M2E, the L2R () and L1R () optimized current flow fields show a factor of about 4.3 and 5.9 times lower current densities, resp. (Table LABEL:Quantification, third column). Overall, L2R and L1R thus have a much higher focality, which can also easily be seen in Figures LABEL:FigDeep-B1,-C1 and -D1. However, if no multi-channel tDCS device is available, the M2E approach provides an optimized bipolar electrode arrangement for a deep and tangential target vector.

The deep target region does not seem to be located in a very deep region of the brain. It gets obvious that the deeper the target vector is located, the higher the averaged maximal current density in the brain compartment, especially in more lateral brain regions. Due to the maximum principle it is thus not possible to target deep regions without stimulating more lateral brain areas. However, many important target regions, e.g. auditory, motor or visual cortex, are located rather laterally, so that it will be possible to target with significant field strength in many applications. Moreover, in many applications of brain stimulation it might also not matter, if non-target regions are also involved, because the experimental setup focuses on the target region, for example when examining the change in event-related potentials (ERP) in pre- and post- tCS stimulation ERP measurements.

## 6 Conclusion and Outlook

A novel optimization approach for safe and well-targeted multi-channel transcranial direct current stimulation has been proposed. Existence of at least one minimizer
has been proven for the proposed optimization methods. For discretization of the respective minimization problems the finite element method was employed and the existence of at least one minimizer to the discretized optimization problems have been shown. For numerical solution of the corresponding discretized problem we employed the alternating direction method of multipliers. A highly-realistic six-compartment head model with white matter anisotropy was generated and optimized current density distributions were calculated and evaluated for a mainly tangential and a mainly radial target vector at superficial locations, an extended target area and a deeper mainly tangential target vector. The numerical results revealed that, while all approaches fulfilled the patient safety constraint, the optimized current flow fields show significantly higher focality and, with the exception of the L2R for the deep target, higher directional agreement to the target vector in comparison to standard bipolar electrode montages. The higher directional agreement is especially distinct for the radial target vector. In all test cases, because of a more widespread distribution of injected and extracted surface currents, the L2R optimization procedure led to relatively weak current densities in the brain compartment. The L1R optimized current density distribution along the target direction was in all test cases stronger than the L2R one and might thus be able to induce more significant stimulation effects. The stimulation will thus enhance cortical excitability especially in the target regions, while it will as good as possible prevent too strong excitability changes in non-target regions.

We were able to demonstrate that the M2E approach provides optimized bipolar electrode montages as long as the target is mainly tangentially oriented. For radial targets, the M2E approach was unsatisfactory, an optimal bipolar electrode configuration might then consist of a small electrode placed directly above the target region with a distant return electrode or a small electrode over the target encircled by a ring return electrode, as proposed in [10].

A further application for the optimization method is transcranial magnetic stimulation (TMS). TMS uses externally generated magnetic fields to induce electrical currents to the underlying
brain tissue [16]. Because there is no safety limit for the total currents applied to the stimulating coils but a safety-threshold for painful muscle twitching [16], the
constrained optimization problem for multi-coil TMS is given as

subject to | ||||

with being the time-dependent magnetic vector potential and = being the threshold for painful muscle twitching [16]. By designing the changes in the magnetic vector potential one can consider as the optimization variables, respectively some parameters on which it depends linearly. The existence of at least one minimizer to the constrained optimization problem for TMS directly follows with similar arguments as in Theorem LABEL:minQE.

Because the optimization method can be applied for both brain stimulation modalities, a combined tDCS and TMS optimization might outperform single modality tDCS or TMS optimizations, similar to what was shown for electro- (EEG) and magnetoencephalography (MEG) [3, 4]. While tCS is able to stimulate a radially oriented target, TMS is mainly not (like MEG is hardly able to detect radial sources [3, 4]). Possible applications of combined tDCS and TMS multi-channel and multi-coil optimization might thus be an improved stimulation of target regions containing both radial and tangential orientations or of deeper target regions. In order to induce action potentials in deeper target regions, the induced current density should exceed the threshold for neuronal depolarization of [16]. On the other hand, the threshold for painful muscle twitching of must be kept [16]. The combination of the optimized tDCS and TMS current density fields might lead to higher current densities in the target and simultaneously reduced current density amplitudes in non-target regions.

While a thorough mathematical analysis of our novel multi-array tDCS optimization method was derived and results for different target regions were presented, besides the first promising results presented in [20], we did not yet further compare our method to the existing approaches in the literature such as, e.g., [10, 30, 21, 28]. Such a comparison is one of our future research goals.

Acknowledgements SW and CHW were supported by the priority program SPP1665 of the German Research Foundation, project WO1425/5-1. MB acknowledges support by ERC via Grant EU FP 7 - ERC Consolidator Grant 615216 LifeInverse and by the German Science Foundation DFG via EXC 1003 Cells in Motion Cluster of Excellence, Münster, Germany.

## References

- [1] B. Agsten Comparing the complete and the point electrode model for combining tCS and EEG, Master thesis, University of MÃ¼nster, 2015
- [2] B. Agsten, S. Wagner, S. Pursiainen and C.H. Wolters Advanced boundary electrode modeling for tES and parallel tES/EEG, submitted for publication.
- [3] Ü. Aydin, J. Vorwerk, P. Küpper, M. Heers, H. Kugel, A. Galka, L. Hamid, J. Wellmer, C. Kellinghaus, J. S. Rampp and C. H. Wolters Combining EEG and MEG for the reconstruction of epileptic activity using a calibrated realistic volume conductor model, PLoS ONE, 9(3) (2015), e93154, doi:10.1371/journal.pone.0093154
- [4] Ü. Aydin, J. Vorwerk, M. Dümpelmann, P. Küpper, H. Kugel, M. Heers, J. Wellmer, C. Kellinghaus, J. Haueisen, S. Rampp, H. Stefan and C. H. Wolters Combined EEG/MEG Can Outperform Single Modality EEG or MEG Source Reconstruction in Presurgical Epilepsy Diagnosis, PLoS ONE, 10(3) (2015), e0118753
- [5] L. J. Bindman, O. C. Lippold and J. W. Redfearn Long-lasting changes in the level of the electrical activity of the cerebral cortex produced by polarizing currents Nature, 10(196) (1962), p. 584–585
- [6] P. S. Boggio, F. Bermpohl, A. O. Vergara, A. L. Muniz, F. H. Nahas, P. B. Leme, S. P. Rigonatti and F. Fregni Go-no-go task performance improvement after anodal transcranial DC stimulation of the left dorsolateral prefrontal cortex in major depression, J Affect Disord, 101 (2007), p. 91–98
- [7] S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends in Machine Learning, 3(1) (2011), p. 1–122
- [8] O. D. Creutzfeldt, G. H. Fromm and H. Kapp Influence of transcranial d-c currents on cortical neuronal activity Exp Neurol, 5 (1962), pp.436–452
- [9] M. Dannhauer, D. Brooks, D. Tucker and R. MacLeod A pipeline for the Simulation of Transcranial Direct Current Stimulation for Realistic Human Head Models using SCIRun/BioMesh3D, Conf Proc IEEE Eng Med Biol Soc. 2012; 10.1109/EMBC.2012.6347236 (2012)
- [10] J. P. Dmochowski, A. Datta, M. Bikson, Y. Su and L. C. Parra Optimized multi-electrode stimulation increases focality and intensity at target, J Neural Eng., 8 (2011), p. 1–16
- [11] J. P. Dmochowski, A. Datta, Y. Huang, J. D. Richardson, M. Bikson, J. Fridriksson and L. C. Parra Targeted transcranial direct current stimulation for rehabilitation after stroke NeuroImage, 75 (2013), p. 12–19
- [12] Eichelbaum, S., Dannhauer, M., Hlawitschka, M ., Brooks, D., KnÃ¶sche, T.R., Scheuermann, G. Visualizing simulated electrical fields from electroencephalography and transcranial electric brain stimulation: A comparative evaluation NeuroImage, 101 (2014) p. 513-530
- [13] R. Ferrucci, F. Mameli, I. Guidi, S. Mrakic-Sposta, M. Vergari, S. Marceglia, F. Cogiamanian, S. Barbiere, E. Scarpini and A. Priori Transcranial direct current stimulation improves recognition memory in Alzheimer disease Neurology, 71 (2007), p.493–498
- [14] F. Fregni, S. Thome-Souza, M. A. Nitsche, S. D. Freedman, K. D. Valente and A. Pascual-Leone A controlled clinical trial of cathodal DC polarization in patients with refractory epilepsy Epilepsia, 47 (2006), p. 335–342
- [15] D. Gilbarg and N. S. Trudinger Elliptic partial differential equations of second order Vol. 224 Springer (2001).
- [16] L. Gomez, F. Cajko, L. Hernandez-Garcia, A. Grbic and E. Michielsson Numerical analysis and design of single-source multicoil TMS for deep and focused brain stimulation, IEEE Trans Biomed Eng, 60(10) (2013), p. 2771–2782
- [17] R. F. Hartl, S. P. Sethi and R. G. Vickson A survey of the maximum principles for optimal control problems with state constraints, Siam Review, 37(2) (1995), p. 181–218
- [18] J. Haueisen, D.S. Tuch, C. Ramon, P.H. Schimpf, V. J. Wedeen, J. S. George, J. W. Bellveau The influence of brain tissue anisotropy on human EEG and MEG NeuroImage 15 (1), pp. 159â166
- [19] R. Herzog, G. Stadler and G. Wachsmuth Directional sparsity in optimal control of partial differential equations, Siam J. Control O Ptim, 50(2) (2012), p. 943–963
- [20] Homölle, S. Comparison of optimization approaches in high-definition transcranial current stimulation in the mammalian brain Master thesis, University of MÃ¼nster, 2016.
- [21] C. H. Im, H. H. Jung, J. D. Choi, S. Y. Lee and K. Y. Jung Determination of optimal electrode positions for transcranial direct current stimulation (tDCS) Phys Med Biol, 53 (2008), p. N219–N225
- [22] J. L. Lions Optimal control of systems governed by partial differential equations Dunod and Gauthier-Villars, Paris (1968)
- [23] M. A. Nitsche and W Paulus Excitability changes induced in the human motor cortex by weak transcranial direct current stimulation, J Physiol, 527 (2000), pp. 633–639
- [24] M.A. Nitsche, D. Liebetanz, A. Antal, N. Lang, F. Tergau and W. PaulusModulation of cortical excitability by weak direct current stimulationâtechnical, safety and functional aspects, Suppl Clin Neurophysiol, 56(3), pp. 255-276
- [25] C. Ortner, W. Wollner A priori error estimates for optimal control problems with pointwise constraints on the gradient of the state, Numer. Math. 118 (2011), 587-600.
- [26] S. Pursiainen, F. Lucka and C.H. Wolters Complete electrode model in EEG: Relationship and differences to the point electrode model, Phys.Med.Biol., 57, pp.999-1017, (2012).
- [27] C. Ramon, P. Schimpf, J. Haueisen, M. Holmes, A. IshimaruRole of soft bone, CSF and gray matter in EEG simulations, Brain Topogr, 16, pp. 245â248
- [28] G. Ruffini, M. D. Fox, O. Ripolles, P. C. Miranda, A. Pascual-Leone Optimization of multifocal transcranial current stimulation for weighted cortical pattern targeting from realistic modeling of electric fields, NeuroImage 89, pp. 216-225, (2014).
- [29] A. Schiela and W. Wollner Barrier methods for optimal control problems with convex nonlinear gradient state constraints, Siam J. Optim., 21(1) (2011), p. 269–286
- [30] R. J. Sadleir, T. D. Vannorsdall, D. J. Schretlen and B. Gordon Target optimization in transcranial direct current stimulation Front Psychiatry, 3 (2012), p. 90
- [31] D. S. Tuch, V. J. Wedeen, A. M. Dale, J. S. George and J. W. Belliveau Conductivity tensor mapping of the human brain using diffusion tensor MRI, Proc Natl Acad Sci U.S.A., (1998). p. 11697–11701
- [32] S. Wagner, S. M. Rampersad, U. Aydin, J. Vorwerk, T. F. Oostendorp, T. Neuling, C. S. Herrmann, D. F. Stegeman and C. H. Wolters Investigation of tDCS volume conduction effects in a highly realistic head model, J. Neural Eng., 11 (2014), p. 016002
- [33] W. Wollner Optimal control of elliptic equations with pointwise constraints on the gradient of the state in nonsmooth polygonal domains, SIAM J. Control Optim. 50 (2012), 2117-2129.