Multi-Region Brain Stimulation Optimization Using Transcranial Direct Current Stimulation

Transcranial direct current stimulation (tDCS) is a type of noninvasive transcranial electrical brain stimulation. By optimizing the current distribution of each electrode on the scalp, the stimulation can be guided to a target brain region using a tDCS dense electrode array system. However, previous studies have yielded simple results using optimization schemes in single target stimulation cases. The detailed parameter settings for each optimization scheme and the associated simulation results have not been comprehensively assessed. In this study, we investigated parameter settings of optimization schemes in detail in both single target and multi-target cases. Two optimization schemes, minimum least squares (MLS) and maximum electrical field strength (ME), were examined in this study. MLS minimizes the squared errors between the expected electrical field and the estimated electrical field, whereas ME maximizes the electrical field strength in the target region. We constructed a five layer finite-element head model with 64 electrodes placed on the scalp according to the EEG 10/10 system for simulation. We evaluated the effects of stimulation using these two schemes under three conditions, 1) single target stimulation, 2) multi-target stimulation, and 3) multi-target stimulation under specific task activation, which shown that directly using MLS and ME scheme in multi-target stimulation case may lead to a wrong result. We also reported the improved results fixed by our proposed weighted MLS and weighted ME schemes which take detailed parameter settings into consideration. Our results indicate that the parameter settings in each optimization scheme greatly affected the final stimulation results, especially in the case of multi-target stimulation, and thus, indicate that the parameter settings of each optimization scheme should be carefully considered according to the expected stimulation mode. Our results also suggest that, by calculating the parameters through our proposed methods, the weighted ME and weighted MLS scheme can precisely distribute energy into each target brain region.

Transcranial direct current stimulation (tDCS) is a noninvasive brain transcranial 2 electrical stimulation method that has been widely used for many clinical and scientific 3 applications, such as epilepsy [1,2], Parkinson's disease [3,4], Alzheimer's disease [5,6], 4 depression [7,8], schizophrenia [9,10] and neuroscience research [11]. tDCS can be used 5 to induce long-lasting alterations of cortical excitability and activity by continuously 6 applying a constant low intensity current through scalp electrodes for a sufficient 7 duration of time [12,13]. Unlike transcranial magnetic stimulation, which generates 8 supra-threshold activation of brain neurons via short-lasting high intensity 9 electromagnetic currents, tDCS modulates spontaneous neural firing in the brain via 10 sub-threshold alterations of membrane potentials [12,13]. Thus, tDCS is a purely 11 neuro-modulatory intervention [13]. 12 In previous tDCS studies, a 0.5-2 mA tDCS stimulation is often applied directly to 13 the scalp via two 20-35 cm 2 electrode pads (one anode and one cathode). The anode is 14 placed on the scalp above the target region and the cathode is often placed on the scalp 15 above the contralateral region, neck, or another area [11][12][13][14]. This stimulation method 16 affects a large brain area between the anode and cathode [14,15]. In some cases, the 17 target region may not even be included in the actual stimulated region depending on 18 the locations of the electrodes. To overcome the drawbacks of this kind of large 19 electrode pads, numerous new tDCS electrode montages have been developed to 20 improve the guidance of the stimulation to the target region. These include 21 high-definition tDCS [16], and dense electrode array tDCS [17][18][19][20][21]. In high-definition 22 tDCS, the stimulation precision is improved by using different types of electrodes or 23 electrode placement. Saturnino et.al. investigated the effects of several types of 24 electrodes and electrode placement on the brain stimulation results, and recommended 25 that using a 4 × 1 ring electrodes placement or two focal center surround ring electrodes 26 can achieve better focality of stimulation [16]. In dense electrode array tDCS, N 27 electrodes or N pairs of electrodes are positioned on the surface of the whole scalp 28 according to a specific principle (e.g., Electroencephalogram (EEG) 10/20 or 10/10 29 system), and the stimulated target is located by using different current distributions on 30 these electrodes (i.e., giving each electrode a specific current). Dmochowski and Guler 31 pointed that when the current of each electrode is calculated using optimization 32 algorithms under certain safety constraints, the stimulation via dense electrode array 33 tDCS can be generally limited to the specific target region [17,18]. 34 When using optimization algorithms, the selection of an optimization scheme (i.e., maximum electrical field strength (ME) inside the target region [18]. In comparison 45 with the MLS scheme, the ME scheme can induce the maximum stimulation strength in 46 the target region and only needs to specify the direction of the expected stimulation. 47 However, these studies only examined optimization performance in the case of a single 48 target stimulation. The performance of the MLS and ME schemes in multi-target 49 stimulation scenarios is unknown. As individual brain functions often involve many 50 brain regions, simultaneously stimulating two or more brain regions may have more 51 impact on a specific brain function. Dmochowski proposed a MLS based optimization 52 scheme that uses the EEG recording to give neural sources which generate this EEG 53 recording a stimulation [20]. Although, this method can give neural sources a relatively 54 precise stimulation without given any additional parameter such as the strength and 55 direction of the expected stimulation, it can only apply to the situation that neural 56 sources which generate the recorded EEG need to be stimulated. Ruffini and his 57 colleagues first proposed the concept of multifocal transcranial current stimulation [21]. 58 In this concept, all brain voxels are stimulated under a specific stimulation plan 59 according to a given pattern from neuroimaging studies, such as a T-map in functional 60 magnetic resonance imaging or a apparent diffusion coefficient map in diffusion tensor 61 imaging [21]. However, this study only proposed a framework. The selection of an 62 optimization scheme, specific parameter settings of each optimization scheme, and their 63 effects on the optimization results were not given in detail.

64
Thus, in this study, we systematically investigated the optimization performance of 65 the MLS and ME optimization schemes when applied to brain stimulation via a 66 64-electrode dense array tDCS system, described previously [17]. We examined the 67 performance of the two schemes in three conditions: 1) single target stimulation 68 optimization, 2) multi-target stimulation optimization, and 3) stimulation optimization 69 under an activation by a specific task.

MATERIALS AND METHODS
Where E represents the electrical field distribution of the brain, V i represents the 80 stimulation potential of the i th electrode area, and N is the number of electrodes. Eq 1 81 indicates that, given each a potential V for each electrode, the potential anywhere in 82 the brain can be calculated by solving Laplace's equation.
T represents the potential of the nodes that are the vertices 92 of these elements, M is the total number of nodes, and K is the coefficient matrix of the 93 linear homogeneous equation system. The K is derived from the space coordinates and 94 corresponding conductivity of each node. Eq 2 not only simplifies the problem described 95 by Eq 1, but also markedly improves the speed at which solutions can be computed, 96 making FEM a popular simulation method in many fields. Solving the potential of each 97 node using Eq 2 enables the calculation of V, the corresponding electrical field E, and 98 the current density J of each node in brain [22].

99
Consider that N electrodes are applied onto the surface of scalp layer of a head stimulation with an arbitrary current magnitude is applied to each anode, the electrical 107 field distribution can be described as follows [17]: Eq 3 to Eq 6 imply that by calculating the distribution coefficient matrix A using FEM, 112 the electrical field distribution E can be calculated as a straightforward linear solution. 113 Thus, given a specific electrical field distribution E, we can use an optimization scheme 114 to calculate a set of stimulation currents I that can generate an E that will maximally 115 approximate E.

117
In this study, we investigated the performance of tDCS in a multi-object stimulation 118 scenario using two optimization schemes: MLS and ME scheme.

119
Minimum Least Squares Scheme 120 The MLS optimization scheme is widely used in many research fields. Given a specific 121 electrical field distribution E, the MLS optimization scheme is defined as follows [17]: Where || · || 2 represents the 2-norm of the vector, e 0 represents the expected electrical follows [22] to minimize the squared error under some specified safety constraints: September 12, 2019 5/24 As decribed in [17], if current satisfies this sum form of current constrains, then 132 |I i | ≤ I max , and | I i | ≤ I max

133
In FEM, the number of nodes M is usually very big (i.e., 10 6 -10 7 or greater). Thus, 134 direct optimization using Eq 8 can be both time and source consuming. Further 135 simplification of Eq 8 is need. The least squares can be extended as: and E T E is a constant. If we let W aa = A T A, w ea = E T A, and w ee = E T E, Eq 8 can 138 be simplified as: As A is already calculated using FEM, E is given, and W aa , w ea , and w ee can be 140 pre-calculated, optimization using Eq 10 can be completed in seconds.

142
In study [18], Guler stated that a MLS scheme should specify the strength and direction 143 of the expected electrical field E. In general, researchers tend to select the direction as 144 the normal or tangential direction [17,18,20]. However, the strength of the expect 145 electrical field e 0 is very important, as it will decide not only the stimulation strength of 146 E, but also the form of E as well (Fig 6). Thus, Guler proposed a new optimization 147 scheme, the maximum current density J inside the ROI, which only requires 148 specification of the direction of the expected electrical field E. As J = σE, where σ is 149 the conductivity of medium, the maximum current density inside the ROI is equivalent 150 to the maximum electrical field inside the ROI, i.e., ME The math formulation of the 151 ME scheme is defined as follows [18]: Where Brain-ROI ||σ(r)E(r)|| 2 2 dr represents the total power outside the ROI that, with 153 regards to safety considerations, should not be larger than p max . · represents the dot 154 product operation. Eq 11 shows that the goal of the ME scheme is to maximize the sum 155 of the electrical field in the target ROI under the safety constraint. According to 156 September 12, 2019 6/24 Guler [18], Eq 11 can also be simplified. The ROI (E(r) d(r))dr can be discretized as: 157 Where V m represents the volume of the mth elements. If we assume uniform 160 distribution of the nodes that construct elements in the brain, then V m can be 161 approximately treated as a constant. Thus, Eq 13 can be approximately represented as: 162 Similarly, Brain-ROI ||σ(r)E(r)|| 2 2 dr can be extended as: Eq 11 can be simplified as [18]: Also, W and Q can be pre-calculated. Optimization operations using Eq 18 can be also 166 completed in seconds.

167
Eq 10 and Eq 11 give the basic formulation of MLS and ME optimization scheme in 168 single target case. However, when we directly used these two schemes in multi-target (i.e. taking all target stimulation brain regions as a whole region), the simulation results 170 were not be the one we desired (Fig 5 & 6). To correctly guide stimulations to target 171 regions in the multi-targets case, we propose the weighted ME and weighted MLS 172 optimization scheme in next two subsections.

173
The Weighted ME Scheme 174 Using the ME scheme described in Eq 18, the W can be rewritten as, Thus, the basic weighted ME scheme can be defined as follows: Where α j represents the weight assigned to the j th ROI and satisfies α j > 0, To account for variations in the sizes of the ROIs, before applying the 178 expected weights to each ROI, the original weight of each ROI, which is mainly 179 proportional to the size of the ROI, should be equalized. Thus, the α j can be defined as 180 follows: Where 1/β j is the equalization coefficient of the j th ROI and S j should be calculated 182 according to specific requirements. Because the size of a ROI is not the only factor that 183 affects the original weight, for accurate calculation, β j can be redefined as: From Eq 22, for accurate calculation of β, single target optimization for each ROI using 185 the ME scheme should be computed first. As this can be time consuming, β as defined 186 in Eq 21 is sufficient for general use.

187
If the goal is to stimulate all ROIs with a similar strength, a strength distribution Where std(·) represents the standard deviation of a vector and s 0 > 0 is the maximum 190 standard deviation. The smaller the s 0 , the more similar the stimulation strength 191 among the ROIs. The main effect of Eq ?? is to constrain the differences between ROI 192 stimulation strengths to a relatively small range.

193
If primary electrodes that maintain the stimulation for each ROI need to be injected 194 with a similar current, the weighted ME scheme should be rewritten as: strength, an additional current distribution constraint should be added: and Where the S j needs to be calculated according to specific distribution requirements, and 203 in each ROI class (e.g., high stimulation class and low stimulation class), S j should be 204 From the MLS scheme described in Eq 10, E can be rewritten as: Thus, the weighted MLS scheme can be defined as follows: September 12, 2019 10/24 and white matter (WM). Briefly, brain tissues were first segmented into scalp, skull, 243 and brain areas using the brain extraction tool (BET) function. Then, the brain area 244 section was further segmented into CSF, GM, and WM using the FMRIB's automated 245 segmentation tool (FAST) function.

246
After segmentation, the five tissues masks were imported into simpleware software 247 (https://www.synopsys.com/ simpleware.html) to generate the final head model (Fig 1). 248 We performed manual correction of the tissue masks using the ScanIP module. Then, this study are detailed in Table 1. stimuli, we predicted steady-state fields without much concern for the 270 charging/discharging effect of tissue capacitance (i.e., using quasistatic solution). The  The ACC locates deeper in brain than PreCG. As shown in Fig 3, for a single target 279 case, the MLS scheme had higher focality (i.e., more concentrated energy distribution in 280 stimulation region) but lower intensity (i.e., lighter color in stimulation region) 281 compared with the ME scheme. The detailed performance characteristics are given in 282 Table 2. From Table 2, we can see that, when the depth of the stimulated region   values gradually approximated those in the ME scheme when power constraints existed 292 outside the ROI (P Brain-ROI ≤ 1e-6 V 2 ·m). For the ROI on the grey matter on the 293 surface of the brain (e.g., PreCG), the performance of the MLS scheme was nearly equal 294 that of the ME scheme when the strength of the expected electrical field was set at the 295 maximum value (5.5 V/m, Fig 4E). The maximum value represents the maximum 296 strength to which the expected electrical field e 0 can be set. Values larger than this will 297 cause optimization results which do not satisfy the power constraint (i.e., 298 P Brain-ROI > p max ). For deeper ROIs in the brain (e.g., ACC), the MLS scheme had 299 better focality but lower intensity compared with the ME scheme when the strength of 300 the expected electrical field was set at the maximum value (1.7 V/m, Fig 4E).

3) When 301
the expected electrical field strength gradually increased to the maximum allowed value, 302 Person's correlation coefficient between the currents calculated by MLS and ME 303 schemes gradually increased to one (Fig 4C). These findings indicate that the     the fifth column in Fig 7). Also, the α in Fig.6 were normalized.  ROI was more concentrated than the one in weighted ME scheme, but the total 347 intensity level was lower than the weighted ME scheme (colorbar in figures). Thus, 348 similar with single target case, the weighted MLS scheme had better focality but a lower 349 intensity compared with the weighted ME scheme.

351
We used activation T-map generated from functional MRI data during a specific space 352 working memory task (SWMT) in 36 healthy humans as a reference for generating a 353 stimulation mode according to the activation strength. The SWMT activation T-map is 354 displayed in Fig 8 and the detailed activation ROI information is given in Table 4. We optimized two stimulation modes using weighted ME and weighted MLS scheme, 356 respectively. For the weighted ME scheme, we used a regional stimulation mode: each 357 ROI had a stimulation strength sum that was proportional to the corresponding peak T 358 value. The detailed parameters for the weighted ME scheme are shown in Fig 9 and the 359 results are shown in Fig 9 and Table 5. The derivation of parameters used in the The optimization results of a regional stimulation mode using weighted ME optimization scheme. The second stimulation mode directly used T values in the SWMT T-map as the 365 weight of each voxel's expected stimulation strength. We kept T values larger than 4.5 366 (p < 0.05, familywise error corrected). The others were set to zero and then normalized 367 using the maximum T value in the map. This stimulation mode was not suitable for the 368 weighted ME scheme, which was regional based. Before optimization, we performed 369 equalization according to the size of each ROI (e.g., number of nodes in the head 370 model). The detailed parameters are shown in Table 6. The result of the optimization is 371 shown in Fig 10 and Table 7 with e 0 selected as 1.5 V/m.   We systematically investigated the optimization performance of two optimization 374 schemes, MLS and ME, using a 64 electrode dense array tDCS system. We first The MLS scheme had been widely used in many brain stimulation studies [17,20]. As 384 its definition can be decomposed into a simple form Eq 10, the processing time can be 385 very short. Dmochowski [17] pointed out the trade-off between intensity performance 386 and focality performance in the MLS scheme. From our investigations (Fig 3 & 4), we 387 found that changing the expected electrical field strength e 0 affected this trade-off.

388
Lower e 0 can increase the focality of stimulation while higher e 0 can increase the 389 intensity. If intensity and focality performance are both important, it may be preferable 390 to search the grid for the expected electrical field strength. However, this operation is 391 known to be time consuming.

392
As an improvement to the MLS scheme, the ME scheme could achieve the maximum 393 intensity performance that MLS scheme could be or might be reached (Fig 4) under 394 power constraints, given in Eq 18. The advantage of the ME scheme is that it is only 395 necessary to determine the expected stimulation direction, while for the MLS scheme, 396 the expected stimulation strength e 0 and direction must be determined. However, this 397 advantage prevented the trade-off between intensity and focality from being controlled, 398 leading to the best intensity performance but the worst focality performance when 399 compared with the MLS scheme (Fig 4). Thus, in the single target case, we recommend 400 the ME scheme for the best intensity performance, and for increased accuracy in 401 controlling the intensity and focality performance, we recommend the MLS scheme.

402
From the Eq 8 and Eq 11, it is the definition that makes MLS and ME scheme have 403 different performance in intensity and focality. The optimization problem of MLS 404 scheme is to minimize the mean square error between simulation electrical field and the 405 desired one. Thus, it can achieve the best focality performance. The optimization can achieve the best intensity performance.

408
Multi-Target Case 409 Normal brain function requires cooperation among many brain regions. Thus, 410 simultaneously applying stimulation to two or more brain regions may be more 411 meaningful in both therapeutic and research applications. We found that directly using 412 the MLS scheme or ME scheme (e.g., taking all target regions as a whole) for 413 multi-target optimization caused the 'overfitting' problem, i.e., the optimization 414 algorithm tended to optimize those regions that were larger in size or had a greater 415 impact on the cost function, and ignore other regions (the first column in Fig 5 and the 416 Fig 6). Thus, we used a weighted MLS scheme and weighted ME scheme to overcome 417 this problem. Our results shown that when applying appropriate weights to each ROI, 418 the weighted MLS scheme and weighted ME scheme could implement an arbitrary 419 regional based stimulation mode (the second column to the fifth column in Fig 5 and 420 Fig 7).

421
For the weighted ME scheme, besides applying weights to each ROI, additional the sum of the electrical field strength was very small compared to the values in Table 451 3. Because ROI 2 was much smaller than L SFG in Table 3   be improved, our proposed methods can still precisely distribute the energy into each 467 target brain region according to any given ratio (Table 5). Thus, our methods can 468 become a guideline for those studies or clinical treatments which need to stimulate two 469 or more brain regions simultaneously using tDCS. genetic algorithms to solve the optimization problem. However, we found that the 477 optimization results of these two algorithms are similar with interior-points method, and 478 more time will be cost by these two algorithms thanrunning the interior-points method 479 100 times. Thus, Future work will focus on the improvement of these two optimization 480 schemes and finding more effective algorithms to solve optimization problems.

481
Secondly, as shown in Fig 3,  This study validates previous works in [17] and [18] in single target stimulation case 492 using MLS and ME optimization schemes, then adds upon their work, by investigating 493 the performance these two schemes in multi-target stimulation case. Our findings 494 suggest that directly using MLS or ME scheme in multi-target stimulation case (e.g.

508
We show the derivation of α in Fig 9. According to the specific stimulation parameters, 510 Eq 25 can be written as: