Measurements and models of electric fields in the in vivo human brain during transcranial electric stimulation

1Transcranial electric stimulation (TES) delivers weak electric currents to the scalp with the goal of modulating endogenous brain activity. Stimulation can be applied as constant unidirectional current (tDCS) or biphasic alternating current (tACS). Weak constant direct currents (tDCS, ≤2 mA) can elicit prolonged shifts in cortical excitability. More recently weak alternating currents (tACS) have also been used in an effort to entrain or modulate brain activity. Because of their simplicity, flexibility and safety profile, these techniques have been investigated in over 70 neuropsychiatric conditions, including major depression, epilepsy, tinnitus, Parkinson’s disease, pain control, and stroke rehabilitation among others. In healthy subjects tDCS may benefit declarative memory, working memory, motor learning, verbal fluency, and planning ability.

2Many clinical trials place stimulation electrodes on the scalp based on the assumption that the ‘active’ electrode will stimulate the underlying brain region with a uniform polarity and irrespective of the placement of the ‘return’ electrode. Computational models of current flow have called this simplistic approach into question. These models capture anatomical details obtained with magnetic resonance imaging (MRI) to predict electric fields inside the brain. Models have shown that the polarity of stimulation is necessarily mixed due to cortical folding. Furthermore, nearby electrodes result in maximal electric fields between the two stimulation electrodes, and not underneath the electrodes as commonly assumed. Detailed computational models also suggest that more focused stimulation (1 cm radius) of the cortex may be achieved by using multiple smaller electrodes. Finally, these models predict that cerebrospinal fluid (CSF) in the inter-hemispheric and ventricular spaces could guide current to deep structures of the brain.

3Unfortunately, despite increasing sophistication in these computational models, none of these model predictions have been directly validated to-date. Early validation efforts for simple spherical models used ex vivo recordings, surface recordings or in vivo recordings in simian. Modern models have compared predictions to human scalp surface recordings, but detailed models have not been validated with intracranial recordings in vivo. Furthermore, the models depend heavily on tissue conductivity values, which have only been measured ex vivo in other animal species. These measures differ from in vivo measurements in humans obtained over half a century ago.

4As a result, there are substantial knowledge gaps in determining the electric fields realistically achieved in the brain. This is of particular concern given that the present estimates (of at most 1 V/m) are at the lower end of what is required to affect neuronal activity in animal experiments. There is also no certainty as to whether the spatial distribution of electric field intensities across the brain is sufficiently well predicted to warrant model-based targeting. Here we aim to address these questions with in vivo intracranial recordings in humans by directly measuring field intensities produced by TES at the cortical surface and deeper brain areas. We report results on ten (10) patients undergoing invasive monitoring for epilepsy surgery, with subdural grids, strips, and depth electrodes. In total we recorded from 1380 intracranial electrodes. These recordings are then compared to various detailed computational models, including differential conductivity between skull spongiosa and compacta, and white matter anisotropy. In doing so, we obtain calibrated models that conclusively answer outstanding questions about stimulation magnitudes, spatial distribution, and modeling choices. For instance, we find that maximal electric field magnitudes are approximately 0.8 V/m when using the generally accepted maximum of 2 mA scalp stimulation. These low field intensities are below the levels typically used in animal studies providing a challenge to future studies on the mechanisms of action of TES. To facilitate validation of new modeling approaches we have made all relevant measurement data publicly available (at http://dx.doi.org/10.6080/K0XW4GQ1).

5Transcranial sinusoidal alternating current was applied to ten patients through two electrodes placed medially over the frontal and occipital pole (patients P03–P011 and P014, Figure 1). For one patient (P014) electrodes were placed at three additional locations (Figure 2A). We first confirmed that measured voltages increased linearly with current intensity (Figure 3A; up to saturation levels of the clinical amplifiers). A modest drop in magnitude was observed with increasing frequency (25% at 100 Hz; Figure 3B), which we ascribe to a non-uniform gain across frequencies for the measurement equipment and electrodes as confirmed with recording of voltages in saline. Thus, we observed uniform gain for the tissue indicating purely resistive medium in this frequency range. Repeated measurements with fixed frequency but minutes of separation indicate the level of stability of recordings over time (Figure 3C). The observed drift in voltage magnitudes over time could have a number of sources, most likely patient movement. Regardless, the variability observed here sets the lower limit in precision one can expect to ±0.11 mV for 1 mA stimulation (standard deviation pooled across electrodes and subjects; Figure 3C, values for P03). Subsequent recordings were performed at 1 Hz and with stimulation current of 0.25 mA to 1 mA, limited by the dynamic range of the clinical amplifier and patient sensation.

6Placement of the invasive recording electrodes was determined by clinical considerations. These electrodes do not necessarily capture electric field magnitudes at the locations with maximal stimulation intensities, nor will they be oriented to best capture the electric field in the orientation of the field vectors. To estimate electric field magnitude and orientation across the brain we built for each subject a realistic anatomical model of conductivity in the head (Figure 4) and calibrated the model by matching measurements with the predicted electric fields (see 'Conductivity optimization'). Recorded voltages and their distribution across the cortical surface are shown in Figure 5A/B for the recorded data and the calibrated, realistic head models. Animated 3D rendering of these images (A-C) are provided as supplementary material for all the subjects (rich media, see Figure 5—source data 1). For this subject (P03) the measured voltages are tightly correlated across locations with the predicted electric potential values. The same is true for all subjects resulting in high Pearson correlation coefficients, r = 0.94 ± 0.04 (Here and in the following all summary statistics indicate mean and standard deviation across subjects). To determine the electric fields at each electrode we calculated the difference of the recorded voltages between neighboring electrodes divided by inter-electrode distance (see 'Voltage and projected electric field measurements'). This provides a measure of the local electric field along the direction of a given electrode pair in units of V/m, which we refer to as the ‘projected electric field’. The measured and predicted values of this projected electric field are shown for one subject in Figure 5E. As expected, the correlation of predicted and measured electric fields is lower than for the raw potentials (here r = 0.86, p = 10−14, N = 45), as the calculated field is the difference of two close-by measurements, each with some inherent noise. Similar results are obtained for all ten subjects (r = 0.81 ± 0.09), suggesting that the spatial distribution of electric field magnitudes is well predicted by the models. When collapsing all recordings across subjects (Figure 5F/G) we find correlation between measured and predicted field projections of r = 0.86 (p = 10−118, N = 405) and r = 0.88 (p = 10−54, N = 164) for cortical and depth electrodes respectively. Note also that the measured field projections in depth electrodes are nearly as strong as in cortical electrodes (Figure 5F vs. Figure 5G; standard deviation of measured field projections across electrodes are 0.059 V/m and 0.065 V/m, respectively).

7Now that we have established the accuracy of the model predictions, we can estimate the electric field magnitudes produced throughout the brain. We ask, in particular, what are the maximum intensities achieved at the cortical surface and in deep brain regions? The electrode configuration used in the experiments gives a range of values within and across individuals (Figure 2B). All values reported here and in Figure 2 are scaled to 1 mA of stimulation. As expected for distant electrodes, strongest electric fields are achieved under the stimulating electrodes (e.g. under the occipital electrode, Figure 2C/D). Maximal field intensities for the mid-line configuration (Figure 1) are 0.28 V/m ± 0.06 V/m across the ten subjects (Figure 6A). For an intensity of 2 mA typically used in clinical trials this would result in fields of 0.56 V/m. For other stimulation configurations on Subject P014 with electrodes away from the mid-line (Figure 2A, which is more typical for clinical interventions), maximum intensity on cortex for 1 mA stimulation was 0.25 V/m and 0.10 V/m at the fronto-lateral and lateral occipital locations respectively. We also model other electrode configurations commonly used clinical trials, e.g., M1–SO (Figure 6A). For the three configurations tested here, maximum field intensities in cortical locations was 0.38 V/m ± 0.09 V/m, again for 1 mA. To provide a sense of what intensities are reached in more extended areas (not just at the 1 mm3 voxel with largest value) Figure 6A also reports the 95th percentile (0.14 V/m ± 0.02 V/m), following previous practice.

8To estimate a margin of error on these values, consider the spread of the measured values shown in Figure 5E. For a given predicted electric field intensity, measurements vary by approximately 0.05 V/m. Maximal recorded values at intracranial electrode locations are about 50% of the maximal predicted values across the entire brain (compare Figure 5F with Figure 6A). Assuming that prediction accuracy scales linearly with predicted electric field intensity, this suggests that prediction error for maximal field intensities is in the order of 0.10 V/m.

9The configuration with both stimulation electrodes on the mid-line (Figure 1) is particular in that significant current is shunted through the highly conducting CSF in the inter-hemispheric fissure. As a result, a relatively strong electric field may be produced in deep brain areas (e.g. peri-ventricular white matter and anterior cingulate, see Figure 2D second row). By measuring distance from the center of the brain ([0,0,0] in MNI coordinates) we can chart the dependence of electric field magnitude with depth for each location in the model brain (Figure 6B). The maximal values at a given depth summarized for all subjects and montages indicate that deep brain areas may experience electric fields that are comparable in magnitude to the cortical surface (Figure 6C, e.g. approximately midway between the center of the brain and the cortical surface (normalized distance 0.25) field magnitude were 0.21 ± 0.04 V/m for these 10 subjects with the mid-line montage).

10In this and the following section we use the empirical voltage measurements to assess the validity of various modeling choices commonly used in the literature. First, conductivity values reported in the literature used in existing models appear to give similar estimate of the electric field magnitudes (see values in 'Conductivity optimization'). For example, measured electric field magnitudes in Subject P03 were close to the predicted values, as indicated by the slope s = 0.72 of the best linear fit of the measured versus predicted field projections shown in Figure 7B. This is true for all subjects (for voltage: s = 0.83 ± 0.24; for field: s = 0.68 ± 0.21. Figure 8C/D under ‘literature’). The conductivity values typically used in the literature appear to correctly estimate the amount of current that crosses the skull and the amount shunted through the scalp. We calibrated the models by adjusting conductivity values for each individual model with the goal of minimizing the mean square error between predicted and measured field projections (see 'Conductivity optimization'; example of a best fit for P03 is shown in Figure 5E). The best fitted conductivity values vary across individuals (Figure 8E–G). Generally the fitting procedure showed that the parameters are not tightly bounded by the data. We estimated the uncertainty of the optimal conductivities by analyzing the sensitivity of the fitting criterion to small changes in these parameters (see 'Conductivity optimization' and Figure 8—figure supplement 1).

11The median of the optimal conductivities are 0.03 ± 0.01 S/m for bone, 0.29 ± 0.10 S/m for skin, 0.38 ± 0.05 S/m for white matter and 0.82 S/m for gray matter (Figure 8E–G; ± indicates the median of the estimation uncertainty across subjects). These differ from the literature values (bone: 0.01 S/m; skin: 0.465 S/m; white: 0.126 S/m; gray: 0.276 S/m), but are largely in the same proportions. During optimization CSF was kept constant at 1.65 S/m. Compared to models using literature conductivities, the models with median values across subjects give significantly better accuracy in terms of predicting the electric field distribution (pairwise t-test: t(12) = 2.36, p = 0.04, Figure 8B), but the magnitude is not significantly different (t(12) = 1.51, p = 0.16, Figure 8D). Note that for magnitude, the t-test was performed on the absolute distance of slope to slope of 1.

12Since the introduction of detailed anatomical models at 1 mm3 resolution there have been significant efforts to improve computational current-flow models by including additional tissue compartments (i.e. CSF, bone layers, and white matter). However, there is an ongoing debate as to the relative merits of each of these model refinements (see 'Guidelines for modeling'). To provide some guidelines to future modeling efforts we built models incorporating various refinements proposed in the literature (Figure 9) and compared model predictions to the empirical data in terms of the spatial distribution of projected electric fields (correlation metric reported in Figure 8B). To test whether individual anatomy is important for current-flow modeling, we compared recordings from one subject to predictions from other subjects. Specifically, for each subject (e.g. P03), the recordings were compared to model predictions obtained from other subjects (e.g. P04–P014) at co-registered electrodes locations. This was done for all the 10 subjects, each one comparing with the other 9 (a total of 90 comparisons). Correlation of recordings with these predictions using a different head model is significantly worse as compared to predicting values from the correct individual model (t-test on the difference between correlations: t(89) = 5.60, p = 10−7). This confirms that modeling individual anatomy is important. Many have argued that it is important to incorporate the entire head anatomy down to the neck. Since the MRI of most patients is truncated at the base of the skull, we extended the field of view (FOV) using a standard head that captures the average anatomy of the lower head (see 'General procedure'). This extended FOV gives significantly better predictions as compared to the original FOV, in terms of the electric field magnitudes (t(10) = 6.17, p = 10-4) and electric field distributions (t(10) = 2.63, p = 0.03; Figure 9B, comparing RMcut vs. RM). Note that P04 and P06 were not included in this comparison as the FOV on these two subjects could not be extended (see 'General procedure'). We next tested the importance of CSF by removing the CSF compartment from the intact model – a model capturing the individual anatomy prior to craniotomy and electrode implantation. This did not significantly affect the distribution of estimated electric fields (Figure 9B, comparing correlations for IM vs. IM-CSF, t(12) = 1.78, p = 0.10), but did deteriorate the accuracy in predicting the magnitude of electric fields (comparing slopes, t(12) = 4.39, p = 10−4). Incorporating heterogeneous conductivities for various bone compartments does not appear to provide a reliable improvement on magnitude estimates (RM vs. RM + 3skull: t(12) = 1.72, p = 0.11), and it seems to give worse prediction on the electric field distribution (Figure 9B, RM vs. RM + 3skull: t(12) = 2.35, p = 0.04). For some patients we had diffusion tensor imaging (DTI) data available (P06, P07 and P09). Previous modeling studies have advocated including such data to capture anisotropic conductivities. Figure 9 shows the result of incorporating DTI information to the realistic model (RM+DTI) including different ways of converting DTI information into conductivity (RM+DTI/VN, RM+DTI/VC, RM+DTI/EIT; see 'Model categories' and Table 1). Three data points are not sufficient for a statistical comparison, but there does not appear to be an obvious advantage to using DTI. Note that all model comparisons were performed using the same fixed conductivity values (from the literature) to avoid biasing the comparisons with optimized parameters. Additionally, model comparisons focused on the accuracy of the spatial distributions, as captured by the correlation values, as they are not strongly impacted by the specific choice of conductivities (Figure 8A/B).

13We set out to conclusively answer questions related to the electric field intensities that can be achieved in the human brain in vivo with transcranial electric stimulation. We also tested whether computational current-flow models accurately predict spatial distributions of electric fields, and which modeling choices are most warranted. Our main finding is that the intensities of electric field reach 0.8 V/m when using 2 mA transcranially. This is close to previous predictions using computational models. Peak intensities are achieved underneath the stimulation electrodes, but also in deep midline structures such as the anterior cingulate and the peri-ventricular white matter for the specific configurations tested here. We calibrated models by adjusting skull, scalp and brain conductivity in individual subjects to provide correct electric field magnitudes. We find that individualized models provide predictions of the spatial distribution of currents with an accuracy of r = 0.86 for cortical electrodes and r = 0.88 for depth electrodes when pooling data across all subjects. These models capture individual anatomy for brain, CSF, skull, air cavities and skin at 1 mm3 resolution. Including details such as anisotropic white matter and inhomogeneous bone compartments does not improve prediction performance. However, extending the FOV to include the entire head and neck significantly improves prediction accuracy.

14There is much debate as to whether electric fields in the brain with transcranial electric stimulation at conventional current intensities (≤2 mA) are strong enough to have a physiological effect. Our study provide a definitive answer to the question of what field intensities can be achieved in human in vivo. Maximal stimulation intensity is difficult to measure directly because the recording electrodes are placed based on clinical considerations, and are thus not arranged where maximal intensities are expected. Subdural electrodes allow measurement of the electric field only tangentially to the cortical surface, whereas strongest field magnitudes underneath the stimulating electrodes are oriented orthogonally to the cortical surface. We thus used the validated and calibrated computational model to estimate the maximal electric field that can be achieved anywhere in the brain. Electric field magnitudes in cortex can be as high as 0.4 V/m for 1 mA stimulation current. For typical electrode configurations used in clinical trials maximal field intensity can reach 0.8 V/m when applying 2 mA. More extended areas can reach a value of 0.28 V/m (95th percentile) under 2 mA stimulation. For some electrode montages, current is shunted into deeper areas through highly conductive CSF with maximal values reaching 0.21 V/m. The present data is broadly in agreement with previous modeling estimates with maximal intensities around 1 V/m (e.g., ) and 95th percentile of 0.35 V/m (e.g., ). While physiological effects in vitro have been found for electric fields as low as 0.2 V/m (for entrainment of coherent gamma oscillations, other effects have required electric fields of at least 5 V/m (to entrain up/down state transitions; see for review on AC effects in vitro). Plasticity effects have been demonstrated only at higher field intensities of ∼20 V/m. It is possible that network effects amplify electric field effects and that larger neurons will be more strongly polarized and thus, in vivo effects in human may be stronger than the effects demonstrated for animal in vitro experiments. Future in vivo animal work may shed light on this fundamental question of the required electric field intensities for various physiological effects (e.g., ).

15A recent study has reported the electric field produced by TES as measured by stereotactic EEG electrodes in two epilepsy surgical patients. While the study does not leverage or evaluate accuracy of computational models, it does estimate maximal electric field to be 0.5 V/m with 1 mA TES. This closely approximates the intensity estimated from our current data (0.4 V/m for 1 mA stimulation). These results are in remarkably close agreement if one considers the differing placement of the stimulation electrodes. We made the deliberate choice of placing electrodes away from the surgical craniotomy to minimize current flow through the highly conducting skull defect. Previous modeling work suggests that electric field estimates do not significantly differ from those for normal anatomy when stimulating electrodes are distant from the skull defect. We confirmed this in our models, which show only minimal differences in field estimates when modeling the craniotomy as highly conducting CSF (in practice it is largely fluid-filled) versus entirely resistive (air filled). However, when electrodes are directly placed over the skull opening, as in, one expects larger electric field intensities than normal as current flows directly into the skull. We thus feel that our estimate is closer to what one may expect in normal skull anatomy. also measured voltages for different stimulation frequencies and found a 10% drop in voltage recording between 1–100 Hz. In our equipment we find a 25% drop which we attribute to a non-uniform gain of the equipment (in saline we measured a gain difference of 25% between 1 Hz and 100 Hz).

16Early pioneering work attempted to validate spherical head models with measurements from a human half-skull suspended in a head-shaped receptacle filled with saline. Predictions of the spherical model were also compared to in vivo recordings from a monkey brain, and scalp measurements in human. In we compared an anatomically detailed model of a human head (of 1 mm3 resolution) with voltages recorded on the scalp surface. But extrapolating the accuracy of model predictions from the scalp surface to the brain is problematic given the ambiguity as to the amount of current entering the skull. There are other indirect efforts to validate models by comparing them to the physiological effects, for instance, motor threshold in vivo in simian and human. performed in vivo intracranial recordings from human, but the data were used to better predict EEG surface recordings, not for calibrating models of electric field with transcranial current stimulation. The most recent studies involve similar in vivo intracranial recordings on human subjects and monkeys, but they do not provide a comparison and validation of current-flow models.

17Accuracy of the models was assessed here in terms of the spatial distribution and the overall magnitude. Accuracy in magnitude was measured as the fraction between measured and predicted values (slope of a linear fit). After calibrating the models the overall magnitudes of electric field were predicted with a single set of parameters across all subjects without a bias (s = 0.84 ± 0.31). Accuracy of the spatial distribution was measured as the correlation of measured and predicted values across recording locations. The individualized models at 1 mm3 resolution provide predictions that correlated with measured electric field intensities (r = 0.81 ± 0.09) across the ten subjects (Figure 8B). These numbers include optimization of conductivities to best match the recordings and may thus be overly optimistic. However, this bias is small as correlation is only minimally affected by the specific choice of conductivities. Indeed, when using a single set of values (median) we find similar performance (r = 0.78 ± 0.10). Estimating spatial distribution of field magnitudes is important if one intends to use models to select the best electrode configuration to stimulate a particular brain region. Many clinical trials and cognitive neuroscience experiments with TES aim to target a specific brain region, and many of the inferences that are drawn from such trials are predicated on having an accurate understanding of which areas are maximally stimulated, what polarity the stimulation has in specific sulci, and which areas are not significantly affected with a specific electrode montage.

18Current-flow modeling approaches range from spherical models to more realistic individualized models derived from MRI. Individualized modeling implies the need to consider the anatomy of individual subjects, in particular for patient populations that may have abnormal brain anatomy. Various teams have also worked on automating segmentation and modeling for this purpose. We confirmed on the present data that modeling individualized anatomy indeed is beneficial in correctly predicting electric field distribution across space. Clinical MRI scans are routinely truncated at the base of the skull limiting the FOV that is available for modeling. By extending FOV with a standard head to capture the anatomy down to the neck we demonstrated that capturing current flows through the entire head significantly improved the modeling results. While most studies argue that the CSF and fluid-filled ventricles should be included in models, others have argued that conventional 1 mm3 MRI resolution cannot measure CSF space accurately and thus, it can be omitted. Our results indicate that including CSF into the model is important for correctly predicting the magnitude of achieved electric field in the brain.

19Other model refinements aim to capture the heterogeneous skull layers. argues that modeling the skull as multiple compartments only has a meaningful effect when there is significant volume of cancelous bone along the current path. A complicating factor is that present segmentation algorithms do not automatically extract cancelous bone, making this modeling particularly labor intensive. In our results, manually segmenting a second bone compartment to account for cancelous bone did not significantly improve results. In addition, some have aimed to incorporate anisotropic conductivities using DTI data. We could not evaluate the benefits of DTI modeling statistically, as we had only three subjects with DTI data. But for these three subjects there was no obvious benefit from adding DTI. This is consistent with previous modeling, which shows limited effects of DTI-derived anisotropy on the estimates of electric field in cortex (most of our electrodes were in cortex). Additionally, there have been debates as to how exactly the diffusion anisotropy (as measured with DTI) should be converted into anisotropy of the electrical conductivity. In fact, the most recent in vivo recordings in humans suggest that white-matter electrical conductivity may actually be isotropic, at least when measured at 50 kHz. Here we evaluated the various approaches proposed in the literature, and found no obvious performance difference between them.

20The overall correlation of predicted and measured fields is 0.86 for cortical electrodes, and 0.88 for depth electrodes (Figure 5F/G). This indicates that despite our best efforts, there evidently remain some model inaccuracies. Simplifying tissue to have uniform conductivity is certainly a source of inaccuracy. Improving on this may require improvements in electrical impedance tomography. A source of uncertainty in our experiment is the placement of the stimulating electrodes. We know that small differences in placement significantly affect model predictions, yet the placement could not be conclusively established as the posterior electrode was occluded by bandages and we therefore only have rough anatomical landmarks that guided the placement. Future efforts may want to establish exact electrode placement using 3D localization devices. Other inaccuracies in the modeling may result from errors in segmentation, in particular thin structures such as CSF or temporal bone, which are in some locations below the 1 mm3 resolution of the model. Finally, openings in the eye socket such as optical foramen and orbital fissure are hard to establish and could serve to guide significant current for the present electrode montages. The models we tested here did not incorporate all the details that could conceivably affect current flow, and indeed, model accuracy could hopefully be improved from the currently reported values by adding such details. For instance, fatty tissue has been modeled by, who suggested that electric fields are altered in particular for obese individuals. Other details that may be important are muscle, arteries, sub-compartments of the eyes, and other anatomical details. A general sense for the sensitivity of the model accuracy to errors in segmentation can be obtained by comparing the prediction performance of the intact models and the realistic models (IM vs. RM in Figure 9). Despite multiple refinements to the segmentation in RM, there was no consistent performance improvement across subjects. To allow for further testing of such refinements we have made the MRIs and recorded voltage values freely available to the public (at http://dx.doi.org/10.6080/K0XW4GQ1).

21Another conundrum in TES modeling is the electrical conductivities of the tissues used for modeling. Values commonly used in the literature are in fact the mean values across multiple references. These conductivities are mostly measured at 10 Hz or higher frequencies. To date the data for the human head measured in vivo under direct current (0 Hz) are rather scarce and date back over sixty years. Just very recently new in vivo measures became available largely confirming the literature values for brain tissue, although these were measured at 50 kHz. Whether these literature values provide correct estimates of electric field magnitude in computational models of the human head has never been directly evaluated in vivo. Here we took a phenomenological approach and adjusted the conductivity of scalp, skull and brain to match the recorded fields at 1 Hz. The obtained conductivities are ‘optimal’ in the sense that they best calibrate the resulting electric field intensities inside the brain. Generally, the model fitting appears to be under-constrained as differing parameters gave similar goodness of fit (Figure 8B). We thus constrained the optimization to only three free parameters. The main goal of optimization was to allow some flexibility for the model to more accurately reflect the overall field magnitudes, which were incorrectly predicted by the conductivity values from the literature. The ‘optimal’ conductivities thus primarily serve to calibrate, on individual heads, the overall magnitude estimates. Once calibrated, these models were then used to estimate maximal field intensities in each head. However, we caution against interpreting these ‘optimal’ conductivity values as the true conductivity values of tissue. These best-fit conductivity values are free to account and compensate as best as possible for all sources of modeling simplifications and errors. For instance, the best conductivity obtained for ‘skin’ has to account for the value of many different soft tissues such as muscle, fat and nerves. This is also true for the skull conductivity, which has to reflect in a single number, different bone layers, which will vary in relative abundance across and within individuals. Also, the optimal conductivities may depend on the location of the stimulation electrodes, segmentation errors, and other modeling choices. Thus, ‘optimal’ values differ considerably across subjects. Nonetheless, the median values we report here appear to be a good compromise for predicting field distributions when individual calibration is not possible.

22To provide a sense of how sensitive model predictions are to choices of the conductivity value we computed the estimation accuracy for the best-fit values. The optimization criterion is fairly narrow around these values with estimation margins for skull, scalp, and white matter conductivities of one tenth to one quarter around the optimal values (Figure 8—figure supplement 1A–C). The model also appears to be sensitive to gray matter conductivity with a 10% change in accuracy for a 10% change in conductivity values. It is important to note that the model fitting process is strongly affected by the electric field magnitude, which is approximately inversely proportional to conductivity values. On the other hand, the spatial distribution of the predicted electric field is much less affected by the detailed choice of conductivities. Thus, it may not be surprising that optimal values differ across subjects by as much as one order of magnitude, yet a single median value provides reasonable accuracy in terms of spatial distributions across all subjects (Figure 8B).

23For the purpose of targeting, the most important aspect of modeling is not the overall magnitude of electric fields but rather their spatial distribution. Targeting aims to select the electrode configuration that achieves highest intensity in one location while perhaps minimizing stimulation intensity elsewhere in the brain. We attempted to provide error bars on the estimated electric fields by comparing repeated measures over an interval of approximately 15–30 min. We found relatively little variation in the voltage measurements (±0.11 mV, e.g., Figure 3C), and note that the variations we did observe were likely due to subject movement. Recent evidence and our own observations confirm that voltage measurements depend sensitively on subject position and even blood pulse. Both phenomena are readily explained by the shunting effect of the highly conducting CSF. Small variations with subject position have been demonstrated to affect the thickness of CSF layer and thus modulate the current flow between brain and scalp. Therefore, the variability of electric fields estimated here (with relatively quiescent subjects) likely underestimates the true variability occurring in practice. We conclude that with r = 0.78 ± 0.10 the spatial distribution of electric field predicted by computational models is likely not far from the actual field distributions when taking into account their inherent variability across time, due to subject motion, position, blood pressure, electrode contact, etc. In this sense, the present study suggests that the current high-resolution model is a good foundation in order to select among different possible electrode configurations, and that model predictions with regard to stimulation focality are meaningful.

24This study was performed in epilepsy patients undergoing surgical evaluation with intracranial EEG electrodes at New York University Medical Center (NYUMC). This protocol was approved by the NYUMC Institutional Review Board (IRB) and all patients provided written informed consent as indicated by the IRB. Eligible subjects were approached and consented generally within 24 hr prior to or after surgery (by author AL). We enrolled 12 subjects between December 2013 to June 2016 (P03–P014). Patients had to be over 18 years old, fluent in English, and able to provide informed consent or have an eligible consenting adult readily available. Subjects were excluded for cognitive impairment (IQ70), a space-occupying intracranial lesion, frequent electro-clinical seizures in the 24 hr immediately after surgery (P012), and incomplete data (P013). Therefore, we have included 10 subjects in this paper (P03–P011, P014).

25Voltages were recorded from implanted subdural electrodes (Ad-Tech Medical Instrument, Racine, WI) – the same electrodes used for intracranial EEG recordings. Electrodes were arranged as grid arrays (8 8 contacts), linear strips (1 8 or 12 contacts), or depth electrodes (1 8 or 12 contacts), or some combination thereof. Intracranial EEG signals were referenced to a two-contact subgaleal strip facing towards the skull near the craniotomy site. A similar two-contact strip screwed to the skull was used for the instrument ground. Grid arrays and strips made from stainless steel and embedded in silastic sheets were implanted on the cortex. They have 2.3 mm diameter contacts with 10 mm center-center spacing. Depth electrodes made from platinum-iridium are inserted deep into the brain (around hippocampus). They have 2.4 mm contact length with 5 mm spacing. The decision to implant, electrode targeting, and the duration of invasive monitoring were determined solely on clinical grounds and without reference to this study. Subdural electrodes covered extensive portions of lateral and medial frontal, parietal, occipital, and temporal cortex of the left and/or right hemisphere (see Figure 1). In total we recorded from 1380 intracranial electrodes: P03 (124), P04 (124), P05 (124), P06 (124), P07 (80), P08 (120), P09 (250), P010 (124), P011 (118), P014 (192). Recordings from grid, strip and depth electrode arrays were made using a NicoletOne C64 clinical amplifier (Natus Neurologics, Middleton, WI), bandpass filtered from 0.16–250 Hz and digitized at 512 Hz. The input impedance of the amplifier is above 100 MOhm, and electrode impedance was kept below 75 kOhm during the recordings. Variation in impedance of up to 50 kOhm across electrodes would contribute at most 0.05 percent error in the voltage measurements. Voltage variations due to varying impedance across electrodes are therefore negligible given other sources of variation, predominantly movement leading to altered current paths and transient noise.

26An hour of the patients intracranial EEG recording prior to stimulation was reviewed by an epilepsy physician (AL) to exclude recent seizures. A physician (AL) performed a pre-stimulation clinical assessment (including assessment of the stimulation skin site and neurological examination) and was present at the bedside during the entire procedure to monitor for clinical safety. Simultaneously, the patients intracranial EEG recording was monitored in real-time at the bedside during stimulation for seizures (AL, DF).

27Consented subjects had one 2 × 2 cm rubber electrode placed on the mid-forehead (Fpz), and a second similar electrode placed at the occiput (Oz) below the sterile dressing and distant from the surgical skull defect, and applied with conductive paste (Figure 1). One subject (P014) had electrode montages placed at differing locations (see Figure 2A). Subjects were then covered with a nickel-cadmium sheath to reduce environmental artifact during recording. Additional sources of environmental noise in the hospital setting were minimized when possible.

28The stimulation protocol used the Neuroconn DC Stimulator Plus (NeuroConn, Germany), with a sinusoidal waveform, at variable frequencies and intensities between 0.25 mA and 1 mA, zero-to-peak amplitudes. Stimulation was performed during waking rest (eyes closed) and afternoon sleep. Maximum stimulation intensities were primarily determined by the threshold at which amplifier saturation occurred. Stimulation was immediately stopped in the event of an electrographic seizure. A repeat clinical assessment (including assessment of stimulation skin site and neurological examination) was performed after stimulation.

29Voltage recordings used the same electrodes as the intracranial EEG. Magnitude and sign of induced voltages were measured by fitting a sinusoid to the recordings during sinusoidal stimulation. This fitting is insensitive to any DC voltage bias of individual electrodes. Except for Figure 3B, all other magnitude estimates are based on recordings during 1 Hz stimulation. These magnitudes were then calibrated by dividing with the intensity of the stimulating current to correspond to 1 mA current injection. The output of this post-processing were plotted and manually inspected electrode by electrode. Those electrodes with strong 60 Hz interference and evident clipping (strong harmonics in the harmonic fit) were discarded and were not used for model validation. In total we were able to use 1205 electrodes among a total of 1380 electrodes: P03 (118/124), P04 (112/124), P05 (118/124), P06 (117/124), P07 (78/80), P08 (116/120), P09 (184/250), P010 (101/124), P011 (111/118), P014 (150/192). Mean and standard deviation per electrode were computed whenever more than one recording session was available (all except Subject P04 and the last two stimulation montages for P014).

30The projected electric field, for both recordings and model predictions, was calculated in the direction of adjacent electrode pairs by subtracting voltage values and dividing by their distance resulting in a measure with unit of V/m. Since electric fields are non-uniform, only local estimates of projected fields were measured by taking the difference only between the closest electrodes. Other field directions and averages could be obtained by looking at more distant pairs, at the cost of making systematic errors due to non-uniformity. Thus, for each electrode, a pair was selected as the closest electrode on the same grid/strip that is within a 10 mm radius for cortical electrodes, and 5 mm radius for depth electrodes. Note that with this operation we are only capturing the magnitude in the specific orientation of the electrode pair used, that is, the cosine with the field vector.

31MRIs were acquired on all the patients prior to the surgical implantation of electrodes. Within 24 hr after the surgery, patients underwent a post-operative MRI to confirm the locations of the implanted electrodes as part of their routine surgical evaluation. The model was built for each patient in the voxel space of the post-op image. However, this image usually has artifacts caused by the magnetic susceptibility of the implanted electrodes. Thus, we instead used the pre-op images to make the models. To this end, the pre-op MRI of each subject was registered and resampled into the voxel space of the post-op MRI.

32The computational models were built following our previous work. Briefly, the registered and resliced pre-op image was segmented by the New Segment toolbox in Statistical Parametric Mapping 8 (SPM8, Wellcome Trust Centre for Neuroimaging, London, UK) implemented in Matlab (R2013a, MathWorks, Natick, MA). Segmentation errors such as discontinuities in CSF and noisy voxels were corrected first by a customized Matlab script and then by hand in an interactive segmentation software ScanIP (v4.2, Simpleware Ltd, Exeter, UK). The field of view (FOV) of the clinical MRI scans was extended down to the neck by applying SPM8 registration functions ‘Coregister’ and ‘Normalise’ to a standard head, and then resampling the registered head and pasting it to the MRI of each individual patient. This process uses the soft-tissue as landmark for registration. It extrapolates reasonably well to align skull and soft-tissue in 8 of the 10 patients. For P04 and P06 we had to leave the FOV restricted to the original MRI scan. The 2 × 2 cm stimulation electrodes were manually created and positioned interactively into the model using ScanCAD module within ScanIP as follows. We took photographs of electrode placement in each patient, and placed the electrode in the model based on the positions shown in these pictures. However, this process was limited by the bandages which occluded the stimulating electrodes, in particular on the back of the head. There the electrode was placed halfway over the occipital ridge unless the photographs suggested otherwise. For each patient a finite element model (FEM, ) was generated from the segmentation data by the ScanFE module in ScanIP and then the electric potential distribution was solved by Abaqus (v6.11, SIMULIA, Providence, RI) under 1 A/m current density injected into the frontal stimulation electrode. The resulting electric potential distribution was calibrated to correspond to 1 mA current injection.

33The intracranial electrodes were localized and reconstructed from pre- and post-op MRIs using previously developed procedures at NYUMC, and their coordinates were registered into post-op images to read off the electric potentials in the brain predicted by the model. Specifically, the voltages on the FEM nodes at the brain that are closest to the electrodes were read out, and are ready to be compared against the actual measurements ('Validation criteria').

34To assess how different parameters in the modeling process affect the model predictions, we built six categories of models for each subject. They are listed in order of increasing complexity in Table 1.

35This is the model that is built on the immediate results from the segmentation process using the pre-op MRI, that is, model with normal head anatomy (gray matter (GM), white matter (WM), CSF, skull, scalp and air cavities). Anatomical alterations introduced by the surgery (e.g., craniotomy) are not included in this model.

36This is built upon the intact model, by further modeling all the structures introduced by the surgery. Specifically, the craniotomy was modeled by inserting a cylinder of the same diameter as the burr hole into the skull and then intersecting the skull voxels. These voxels were given CSF conductivity and this area was covered by soft-tissue (skin conductivity). Additional instruments used during the surgery, such as the subgaleal electrodes and the Jackson-Pratt (JP) drain used to vacuum blood products, were modeled by manually drawing out the corresponding voxels from the post-op MRI. Since the intracranial electrode contacts are extremely small (the thickness is less than 1 mm in grid arrays and linear strips), we did not model them explicitly. Instead, we did model the geometry of the grid arrays and linear strips, by taking voxels from the CSF. To this end, customized Matlab script was written to generate cylinders along adjacent intracranial electrodes on the strips and arrays and then intersect with the CSF voxels. The insulated wires running through the scalp were not modeled as their diameter (<1 mm) is not expected to significantly affect the predicted current flow. Figure 4A–F shows a 3D rendering of such a realistic model for Subject P06. The realistic model also includes an extension of the field of view as described above.

37This is the same as the realistic model, except the lower head is cut off due to the limited FOV of the clinical MRI scans. Note for Subjects P04 and P06, RMcut is the same as RM, since for these two subjects we did not extend the FOV when building the models (see 'General procedure').

38To test whether CSF is really important as claimed by many modeling studies, we removed the CSF in the intact models by combining it with the brain (eyeballs were combined with the scalp).

39Starting from the realistic models, we further differentiated the skull spongiosa and compacta. Specifically, manually determined thresholding was performed for each subject to segment out the spongy bone from the compact bone. See Figure 4G for an example. Note that we did not make the skull anisotropic, as three-isotropic-layered model was recommended by a modeling study and ex vivo measurement. Thus, this skull model is inhomogenous and isotropic.

40Anisotropy (conductivity varying with direction) was added to the gray and white matter on three subjects for whom diffusion tensor imaging (DTI) scans were available (P06, P07, P09). The DTI scans were first corrected for non-linear susceptibility and eddy-current distortions with the FSL ‘topup’ and ‘eddy’ command-line tools respectively. The mean of the non-diffusion-weighted (b0) volumes was then used to create the DTI brain mask via the FSL ‘BET’ utility. Since the conductivity tensor is required to be positive semi-definite, the diffusion tensor was calculated for the brain voxels with a non-linear optimization relying on positive-definiteness constraints as implemented in the Camino software. After this, to integrate the DTI data into the model, the fractional anisotropy (FA) of the diffusion tensor was computed and registered with an affine transformation to the voxel space of the post-op MRI (where the model is built). To circumvent the complexity of spatially transforming the diffusion tensors, the previous transform was applied to all the distortion-corrected diffusion-weighted volumes. The positive semi-definite diffusion tensor was then computed again from the transformed diffusion-weighted data, but within a different brain mask defined by the model. See Figure 4H for a slice view of the diffusion tensor of Subject P06. Each mesh node in the model was then associated to the nearest diffusion tensor. These diffusion tensors were then converted into conductivity tensors using one of four different algorithms used in the literature: direct mapping; volume normalized; volume constrained; and equivalent isotropic trace. Models with these anisotropic conductivities were also solved in Abaqus.

41Note that when comparing across different categories of models listed above, literature conductivity values (see 'Conductivity optimization') were used when solving the models. This is to avoid any potential bias that may be introduced when using optimized conductivities ('Conductivity optimization').

42For both the voltages and the projected electric fields, the model-predicted values and the recorded values (see 'Voltage and projected electric field measurements') were compared for each subject in terms of two criteria: (1) Pearson correlation coefficient between recorded and predicted values; (2) the slope of the best linear fit with predicted value as ‘independent’ and measurement as ‘dependent’ variables. The correlation captures how well the model predicts the distribution patterns of the electric potential regardless of a potential mismatch in overall magnitude. The slope measures how accurate the model estimates the absolute magnitudes. Since the ground reference for the recordings (the subgaleal electrodes) is different from the ground for the model (electrode placed on the back of the scalp surface), mean across electrodes was subtracted from the recordings and the model outputs before they were compared. Slope and correlation are both independent of an overall offset.

43Pairwise t-test was carried out to compare how different conductivity values and different model types affect the prediction accuracies (Figure 8 and Figure 9). Note that when comparing magnitude accuracy as indicated by the slope, the t-test was performed on the absolute distance of to slope of 1.

44We first assigned literature conductivities to the realistic models. The following values were used (in S/m): GM–0.276; WM–0.126; CSF–1.65; skull–0.01; scalp–0.465; air–2.5 10−14; gel–0.3; electrode–5.9 × 107. For realistic models with inhomogeneous skull: spongiosa: 0.02865 S/m; compacta: 0.0064 S/m. The JP Drain and intracranial electrode grids and/or strips were made almost insulated (<10−14 S/m).

45We note that in the literature cited above a range of values are reported, all of which have been measured at 10 Hz or higher frequencies (up to 1 GHz, in S/m): skull: 0.0028 ∼ 0.08; scalp: 0.0003 ∼ 1.0; WM: 0.1 ∼ 1.0; GM: 0.016 ∼ 1.0.

46Secondly, the optimal conductivity for each patient was obtained by fitting the model outputs to the recordings. Formally, the following optimization problem was solved:

47Here is the projected electric field on electrode calculated from recorded voltage, and is the corresponding model-predicted value, which depends on the conductivity σ. M is the total number of intracranial electrodes. Mean across electrodes was subtracted from both recordings and model predictions before evaluating the cost function. A general purpose optimization algorithm – Pattern Search, available in Matlab – was used to solve this optimization. This procedure was carried out only on the realistic model for each subject individually. In fitting conductivity parameters we found that the values are generally under-constrained (very different conductivity values can give similar goodness of fit). We therefore limited the fitting to adjusting scalp, skull, and brain conductivities, but fixed CSF (1.65 S/m) and maintained the ratio between gray and white matter constant (=2.19) reflecting the values reported in the literature. Thus there are only three free parameters in the optimization, which is also beneficial in order to keep computing times reasonable. Each model evaluation took approximately 30 min, and optimization was terminated after 300 evaluations of the Pattern Search algorithm. Note that the conductivities of the tissues are optimized jointly, that is, they are adjusted simultaneously per evaluation of the cost function.

48To estimate the uncertainty in this parameter estimation we computed the variance of the estimator based on the Cramér-Rao bound. This was done separately for each conductivity value, which gives a loose lower bound for the estimation error, that is,

49where is the estimation variance of the conductivity, and is the residual error (Equation 1). Here is the parameter of a quadratic fit to the cost function around the optimal value. The expression to the right is the inverse of the Fisher information, estimated numerically. For parameters that were fixed to the literature values (gray matter and CSF), we similarly solved the model for close-by values (±10%) and report the sensitivity of the normalized error to this 10% fluctuation. All values are shown in Figure 8—figure supplement 1.

50Finally, we also tested the median of the optimal conductivity across all subjects (P03–P011, and P014).