Finite Element Modelling Framework for Electroconvulsive Therapy and Other Transcranial Stimulations

Electroconvulsive therapy (ECT) is widely acknowledged as a highly effective treatment for major depressive disorder, and transcranial brain stimulation techniques in general are of great interest for therapeutic neuromodulation and neurostimulation. It is however difficult to determine the effect of electrical stimulation on the brain due to the complex current pathway between the electrodes, which cannot be readily visualized. Computational models of the human head, combined with a finite element implementation of the Laplace equation, can be used to provide information on the electrical stimulus, such as voltage, current density and electric field distributions, helping to understand the effect of transcranial stimulation on particular brain regions of interest. In this chapter, a detailed protocol for creating a finite element computational head model for transcranial electrical stimulation is provided. Procedures outlined include image segmentation, white matter anisotropy extraction, meshing and finite element model implementation. The computational modelling methods described here can be used, for example, for future novel designs of improved ECT protocols.


Introduction
Electroconvulsive therapy (ECT) has been used to ameliorate major depressive disorder for patients who are resistant to drug therapy. The treatment involves applying a train of alternating pulses across two electrodes placed on the scalp. ECT is an effective treatment [1], but also carries a risk of cognitive side effects, such as disorientation and memory loss [2]. Treatment efficacy has been noted to rely on multiple factors, such as electrode placement and stimulus dose [3]. In addition, there is currently also great interest in other brain stimulation techniques for therapeutic neuromodulation or neurostimulation, including transcranial direct current stimulation.
Due to electrical conductivity variation across different tissues, the current pathways induced by electrical stimulation are not straightforward to identify. The presence of highly resistive skull and air-filled paranasal sinuses impedes the passage of electrical currents, forcing the majority of currents to travel through the less resistive regions [4]. Furthermore, the white matter exhibits a strongly anisotropic conductivity due to its myelinated structure, thus determining the prefential current pathway within the brain [5,6]. As such, the electrical current distribution in the brain resulting from brain stimulation is complex and cannot be readily imaged, and is impractical to be measured empirically. Alternatively, the electrical current and electric field (E field) in the brain can be simulated via computational modelling. The finite element (FE) method is one of the most popular numerical approaches for solving models expressed as partial differential and integral equations.
The main goal of computational modelling for ECT and other brain stimulation techniques is to determine the region(s) modulated by the electrical stimulus. It is believed that non-invasive brain stimulation shifts the tissue's membrane potential, subsequently affecting neuronal firing [7]. The use of computational modelling to examine differences in regional E fields as the ECT stimulation approach is altered allows for a better understanding of the relationship between brain stimulation and clinical effects with current forms of ECT, as well as offering the potential for futher improvements in ECT stimulation techniques. In this chapter, we will discuss approaches and steps necessary to implement computational modelling of the human head to determine the voltage and E field distribution during the application of ECT and other transcranial electric stimulation techniques. Figure 2.1 describes the steps needed to undertake a computational study of electrical brain stimulation. Similar steps have been performed as part of previous finite element studies [3,5,8,9], with minimal variation to suit the need for each study.

Image Pre-processing
In order to simulate the properties of different head structures, these structures need to be individually reconstructed from the acquired images. The process of partitioning the image into different domains or "masks" is known as image segmentation.
In order to increase the accuracy and reduce the effort of image segmentation, certain pre-processing procedures are performed prior to segmentation. These may include resampling (to reduce the resolution), cropping (to restrict the image set to the region of interest, i.e. ROI), artefact correction (such as motion, metal and bias field artefacts), edge and contrast enhancement and image registration. These operations can be performed using a selection of open-source image-processing software, such as ImageJ (https://imagej.nih.gov/ij/), 3D Slicer (https://www.slicer.org/) and ITK-SNAP (http://www.itksnap.org/). Among these, bias field correction and image registration are highly common pre-processing steps in the segmentation of MR head scans.

Bias Field Correction
Bias field noises are caused by low intensity and smooth signals that distort the MRI images, and are present especially in older MR devices [10]. This type of noise causes regional differences in signal intensity in the images, leading to non-uniform intensities in the same head structure, as shown in Fig. 2.2. If left uncorrected, segmentation quality may be affected. Bias field correction can be performed prior to segmentation using open-source tools designed specifically for head segmentation, as listed in Table 2.1.

Image Registration
It is not uncommon to obtain multimodal MRI scans, such as T1-weighted scans together with T2-weighted, proton density (PD)-weighted or diffusion-weighted scans, to provide complementary information regarding tissue structures in the brain. As these scans may be acquired in different coordinate systems, it is essential to perform image registration to transform these into the same coordinate system prior to segmentation. A common registration method is affine transformation, which is a linear transformation aligning two sets of images together through translation, scaling, shear mapping and rotation [11]. When a linear registration method is not able to provide a satisfactory outcome, such as when registering scans from different subjects, a non-linear transformation method should be applied. These transformations are performed, automatically or manually guided, through identification of anatomical landmarks, such as the corners of the ventricles, which are easily distinguishable from the images. Image registration is typically available in image-processing software packages.

Image Segmentation
The structural domains are separated through their identifiable landmarks and/or edges. T1-weighted MRI is typically used as it provides a good contrast between whole head structures, especially between grey and white matter. Skull extraction can be challenging, especially at the ethmoid sinus region, but this can be rectified by combining the T1-scan data with CT, T2-weighted or PD-weighted MRI scans [12].
Several open-source software packages designed for brain segmentation have been developed by different research groups. These can automatically extract major head structures, such as grey and white matter, skull and cerebrospinal fluid (CSF), from MRI images. Several of these can also further partition the grey matter into various cortices based on pre-defined atlases. A list of automatic brain segmentation software packages is provided in Table 2.1. It is good practice to perform additional checks following automatic segmentation to ensure there are no segmentation errors. This is usually performed in image processing software that allows manual segmentation, e.g. open-source 3D Slicer and ITK-SNAP, as well as commercially available tools such as Materialise Mimics Innovation Suite (https://www.materialise.com/), Amira (https://www.fei.com/software/amira-for-life-sciences/) and Simpleware (https://www.synopsys.com/simpleware). Other processing, such as smoothing, can also be performed to improve segmentation quality. Several software packages provide training datasets and online tutorials to assist learning.

Manual Segmentation
Thresholding is a critical step in manual segmentation. A thresholding filter can be applied to select particular brain regions. For example, grey and white matter can be easily discerned from T1-weighted MRI images since they exhibit different image intensities. As such, they can be readily segmented into individual masks. In addition to grey and white matter, the CSF space can also be easily recognised from its high intensity in T2-weighted images. This facilitates masking of CSF in between the pial surface (outer grey matter surface) and the dura surface (inner skull surface) as well as the interior brain ventricular system. Segmentation can also be performed with other image processing techniques such as "seeding and growing", Boolean operations and mask growing/shrinking. These options are available as standard features of most image segmentation software packages [13]. "Seeding and growing" begins by manually placing seed points in a particular region. The seed points are then expanded to adjacent pixels based on certain region membership criteria, such as pixel intensity and connectivity, until all the connected pixels cover the structure of interest. This prevents the inclusion of other regions of similar pixel intensities into the same mask: for example, segmenting the brain by thresholding alone, ignoring connectivity, may inadvertently include the bone marrow of the skull.
Boolean operation techniques work directly on segmented masks, creating a union, intersection or difference between two masks. These can be used to obtain regional domains encapsulated between two domains. For example, the CSF is encapsulated between the dura and pia structures. Rather than directly segmenting the CSF, it can instead be obtained by performing Boolean subtraction of the encapsulating domains enclosing the brain and skull. This ensures continuity of the surfaces between domains in addition to segmentation efficiency. Boolean operations can also be used to detect boundary intersections between masks resulting from segmenting errors, as detailed in Section "Challenges and Tips in Segmentation".
Mask growing/shrinking is another technique that operates directly on segmented masks. It is similar to performing a mask scaling. Depending on the algorithm, this may be performed in 2D or 3D, or using a uniform or non-uniform approach based on pixel intensity. The combination of these two operations may be used to remove islands, close holes, or interpolate between every two or three image slices.

Surface Smoothing
Following manual or automatic segmentation, output masks are often rough and contain sharp edges. Small islands, i.e. disconnected shells, may also be formed during segmentation and should be removed. These issues can be addressed using a smoothing process as shown in Fig. 2

.3.
Surface smoothing can be performed on the masks within the image processing software, using either Gaussian, median or Laplace smoothing. It can also be performed after the masks have been exported as surface triangulated objects, usually in .stl (stereolithography) format. Operating platforms that can perform smoothing include Blender (https://www.blender.org/), Geomagic Wrap (https:// www.3dsystems.com/software/geomagic-wrap) and Materialise 3-matic (https:// www.materialise.com/). The smoothing strength must be tuned so that the accuracy of the structure is not compromised. Any sharp edges in the form of spikes need to be removed, since this may prevent efficient meshing in later stages of the modelling effort. Furthermore, such a structure is unlikely to be correct, particularly if located between the brain gyri. For this reason, it is often of interest to observe the effect of ECT or other transcranial electric stimulation techniques on specific brain substructures. Nonetheless, extracting these structures can be challenging since it requires knowledge of anatomical landmarks, which are not readily discernible. In addition, manual segmentation can risk losing consistency among multiple subjects, especially if the segmentation is performed by different people. Some brain segmentation software provide automatic labelling of brain regions, as shown in Fig. 2.4. This labelling is based on established brain atlases where each region has been meticulously mapped. Examples of brain atlases are BrainSuite's BCI-DNI_brain [14] and USCBrain [15], which are available in the BrainSuite software.
The first step in brain region labelling involves registering the subject's skullstripped brain images to the atlas. This utilises linear and nonlinear warping to align the subject's brain with the atlas. Subsequently, the regions are automatically labelled based on anatomical landmarks [14,16].

Challenges and Tips in Segmentation
Several challenges can be encountered during the segmentation process. In this section, several tips and precautions are presented. Unwanted intersections between surfaces can present difficulties during the volumetric meshing stage. As these intersecting surfaces are usually small, as shown in Fig. 2.5, very fine mesh elements will be created around these surfaces. The intersections could also lead to errors during simulation since these structures are merely segmentation artefacts. For example, if the inner skull surface protrudes into the outer skull, this will create a hole in the skull domain, connecting the CSF directly to the scalp. Since the skull is highly resistive, a preferential current pathway is unintentionally created, producing a simulation error.
To prevent such intersections, segmentation can be performed from the outermost surface first, gradually moving towards the inner surface. Inner surface segmentation can then be performed using the outer surface boundary as a guide. Small intersections can also be highlighted using Boolean operators and subsequently corrected manually. It should be noted that intersections can exist even for automated segmentation tools such as those of FSL and BrainSuite. As such, the segmentation output must always be visually inspected.
Other challenges can also arise when segmenting the skull in T1-weighted images, since both compact/cortical bone and CSF appear dark. This may be resolved by directly obtaining the skull mask from CT scans. However, it is often difficult to acquire CT scans due to concerns over unnecessary radiation exposure. The skull is thus extracted by defining the outer skull surface using T1-or PD-weighted images, and the inner skull surface using T2-or PD-weighted images. The spongy bone of the skull is often identifiable as a bright region between the two thin dark regions (compact bone) outside the brain. Segmenting the air-filled paranasal sinuses within the skull may also be difficult as they appear dark. They are often recognised as frontal bone regions where spongy bone is missing.
It is not uncommon to generate a mutually exclusive mask for each domain. This may however create surface continuity problems, especially after surface smooth-

Fig. 2.5
Intersections between grey and white matter due to segmentation error. This can be overcome if segmentation of the grey matter is performed first, followed by white matter segmentation ing, during which both inner and outer surfaces of each mask are exclusively modified. A better practice would be creating inclusive masks, i.e. an outer domain mask within which all other domains of interest are contained. For instance, instead of creating a grey matter-only mask, a brain mask that also contains CSF in the ventricles as well as white matter can be generated. Separation can be accomplished using Boolean operations after all necessary modifications are performed. In some FE meshing tools, separation is not only unnecessary but may also cause contact surface problems between domains.

White Matter Anisotropy
White matter consists of myelinated neuronal tracts, which contribute to its highly anisotropic behaviour. This microanatomical characteristic influences the electrical conductivity such that it is more conductive along the tract than in the transverse direction [6]. Consequently, this affects the spread of the ECT electric field. Simulation studies by Lee et al. [5] and Bai et al. [6] show that disregarding this anisotropy can result in errors in deeper brain structures such as the corpus collosum and hippocampus. As such, the anisotropy helps direct current towards deeper brain regions, where significant effects have been noted following ECT [17].
The linear relationship between electrical conductivity and water diffusion tensor has been experimentally validated [18,19], suggesting that the conductivity tensor shares the same eigenvectors as the diffusion tensor. The water diffusion tensor can be extracted using the diffusion tensor model for diffusion-weighted MRI (DW-MRI) [20], which can be performed in FSL using the probabilistic tracking algorithm from the FDT diffusion toolbox [21][22][23].
Two separate files containing b-values and b-matrices for all gradient directions are required as input for the diffusion tensor calculation. The former summarises the sensitivity to diffusion for each gradient direction, whereas the latter reflects the attenuation effect in x, y and z for each gradient direction [20]. In addition, the input also requires a 3D NIfTI image file of the brain region of interest (ROI), and a 4D NIfTI image file combining all gradient direction scans. Eigenvectors and fractional anisotropy (FA) are then calculated.
FA determines the anisotropic characteristics of a tensor. In brief, an FA value of 0 indicates complete isotropy whilst a value of 1 indicates complete anisotropy. The FA value is determined from the eigenvalues (λ 1 , λ 2 , λ 3 ) of the diffusion tensor as follows [24]: wherel is the average of the three eigenvalues. Regions with a low FA value (typically FA < 0.45), which suggests a low local anisotropy, are removed from the tensor analysis [25] (Fig. 2.6). The conductivity tensor of white matter, σ, is calculated from: where S is the orthogonal matrix of unit eigenvectors obtained from the white matter diffusion tensor, and σ l and σ t are the conductivities in the longitudinal and transverse fibre directions, respectively, which may be calculated using various methods [26,27]. Following the diffusion tensor calculation, conductivity tensors of data points in the DW-MRI scans can then be linked to their individual coordinates. Only conductivity data with a strong anisotropy signal (FA ≥ 0.45) should be exported. Afterwards, the conductivity at the undefined region can be linearly interpolated using the neighbouring strong anisotropy signals.

FE Meshing
After segmentation, the masks are exported as triangulated surface objects, usually in the form of an .stl file. The masks need to be polyhedralised (tetrahedralised in most scenarios), before they are ready to be used in FE analysis. Polyhedralisation (or FE meshing) is the process of generating polyhedral mesh elements to approximate a geometric domain, and these elements are the basis of the FE method. FE meshing is performed in dedicated meshing software, such as the open-source SALOME (https://www.salome-platform.org/), Materialise 3-matics, Simpleware, as well as ICEM CFD and Fluent which are both from ANSYS, Inc. (https://www. ansys.com/). The volumetric mesh will then be imported into FE simulation software such as COMSOL Multiphysics (https://www.comsol.com/) or ANSYS Workbench. A tight contact between masks is essential to ensure continuity between meshed domains. It is thus often advisable in many of these meshing software packages to not import perfectly mutually exclusive masks. One practice is to generate inclusive masks, as detailed in Section "Challenges and Tips in Segmentation". Another is to generate an open intersecting surface in one mask; for example, an open intersecting surface at the end of the spinal cord extending beyond the closed surface of the skull. Afterwards, an intersecting curve must be formed at the intersection so that both surfaces can be meshed, followed by the remaining volumes.
In addition, the shape of these triangulated surface objects may be approximated with non-uniform rational basis spline (NURBS) surfaces, whose shapes are determined from a series of control points. The NURBS-approximated objects, exported in IGES format, can be imported as geometry into FE software and readily meshed into a cluster of polyhedral mesh elements. This provides flexibility in modifying the model geometry directly within the FE analysis platform. The NURBSconversion process is available in 3-matics, Geomagic, Blender etc.; however, tedious manual operations are inevitable if the object has a complex structure.
Whilst meshing is mostly performed automatically by such software, care should still be taken in setting up the meshing method, including maximum and minimum element size, element growth factor, mesh smoothing parameters and if necessary, mesh coarsening paremeters. It is also recommended to check mesh quality to identify poor quality elements, duplicate elements or uncovered faces. These errors may need to be repaired manually.

Physics and Property Settings
Bioelectromagnetism is the study of electric, magnetic and electromagnetic phenomena arising from living cells, tissues or organisms. In the field of bioelectromagnetism, biological tissues are generally considered as "volume conductors", in which the inductive component of the impedance is neglected, and resistances, capacitances and voltage sources are distributed throughout a three-dimensional (3D) region [28].
In the low-frequency band, where the frequency of internal bioelectric events lies, capacitive and electromagnetic propagation can be neglected [29,30], thus treating bioelectric currents and voltages in living tissues as stationary [31]. This is known as the quasi-static approximation. A recent modelling study by Bossetti et al. [32] Fig. 2.6 (a) Original diffusion tensor images with colour denoting the principal direction (largest eigenvalues) of the diffusion tensor (red indicates left-right, green indicates anterior-posterior, and blue is inferior-superior). (b) Fractional anisotropy with lighter colour indicating higher anisotropy. These images were generated using 3D Slicer software investigated the difference in neural activation between solving the quasi-static field approximation and solving the full inhomogeneous Helmholtz equation using square-pulse current stimuli. They found that for commonly used stimulus parameters, the exact solution for the potential (including capacitive tissue effects) can be well approximated by the quasi-static case. Given the relatively low values of permittivity and magnetic permeability in living tissues, the quasi-static approximation can therefore be employed in computational head models of transcranial stimulation.
The electrical potential φ resulting from ECT can be obtained by solving the Laplace equation: where φ is the electric potential and σ is the electrical conductivity tensor. In order to solve Eq. 2.3, boundary conditions have to be defined at all domain boundaries/ surfaces. Based on the "quasi-uniform" assumption, the degree of activation in a target region is proportional to the local electric field magnitude  E [33]. The electric field  E is determined from Maxwell's equations under quasi-static conditions using Simulation results can be analysed by comparing the average electric field magnitude E in several ROIs of the brain, which is determined using: where  E is the local electric field magnitude at every spatial point in the ROI, and the denominator is simply the volume of the ROI.

Tissue Conductivity
It is typical to assume homogeneity and isotropy in most head tissues, except for white matter. Tissue conductivities for each domain, presented in Table 2.2, were determined from previous experimental studies, as described in Bai et al. [25].

Electrode Placement
Conventional ECT electrode placements including bifrontal (BF), bitemporal (BT) and right unilateral (RUL) placements have been substantially investigated using computational modelling [3,5,26]. Several variations of electrode placement have also been investigated with estimates of electric field strengths in key brain regions, with an aim to improve existing ECT protocols [34,35].
There are several ways to define the location of an electrode on the scalp: 1. Creating a geometric object representing the physical ECT electrode on the scalp, e.g. a short cylinder with a diameter of 5 cm and a conductivity of 9.8 × 10 5 S/m [5]. A normal inward current density J| electrode can be defined at the top boundary of one electrode, where I is the stimulus current, e.g. 800 mA for ECT, and A is the top boundary electrode area. A ground condition (φ = 0) can be defined at the top boundary of the other electrode. Current continuity or other conditions representing the skinelectrode interface can be defined at the scalp surface of both electrodes.
2. An isolated geometric boundary (e.g. a circular boundary of diameter 5 cm) defined on the scalp surface. The boundary condition for the ECT electrodes can thus be defined, with other conditions representing the skin-electrode interface over this isolated boundary. 3. A mathematically defined boundary created by intersecting the scalp surface with a geometry defined by analytic functions (e.g. an analytically defined sphere with a diameter of 5 cm for creating ECT electrodes) [25,36]. An evenly distributed normal current density J is then applied over the analytically defined geometry. A normal inward current boundary condition is defined over the entire scalp such that everywhere, except at the ECG electrodes, the normal current density is zero. The other electrode is defined to have a normal outward current density, −J| electrode . Note that setting the conductivity to zero in any tissue region may present a numerical issue in the simulations. The common practice is to either set this domain as inactive in the simulations or to assign an extremely low value such as 1e-8 S/m

Other Boundary Conditions
The neck region on the lower boundary of the head can be defined using a ground [25] or distributed impedance boundary condition [3]. The latter sets the neck as being connected to a resistive compartment, which in turn is connected to ground. This permits some flow of current through the neck boundary, albeit negligible in magnitude (0.001% of delivered ECT current). The boundary condition of this distributed impedance is defined as follows: where  n J × neck is the normal outward current density at the neck boundary, σ neck is the conductivity assigned to the boundary, d s is the thickness of the boundary and φ ref is the reference voltage, which is set to zero for ground.
The rest of the scalp is set as an insulated boundary (i.e. zero normal component of current density). If the sinuses are inactive in the simulation, their boundaries should be set to insulated as well. All other internal boundaries must be set as continuous current density interfaces to ensure electric continuity between domains.

Numerical Solver Settings
Under quasi-static assumptions, the stimulus amplitude and the resulting voltage, electric field and current density are all in a linear relationship. Therefore, it is sufficient to employ a steady-state solver. A detailed head model usually involves computing over a large number of mesh elements (>5 million). In COMSOL, there are two classes of linear solvers for computation: direct solvers, such as PARDISO, which are time-efficient but require large computational memory, and iterative solvers, such as conjugate gradient, which approach the solution gradually, and thus are memory-efficient but may require substantial computational time. An iterative solver is a better choice for standard desktop workstations with 24 to 64 Gb RAM. The absolute tolerance of the error in previous works was set to 10 −5 [36] or even at 10 −8 [5,9]. In general, a lower absolute tolerance yields a more accurate result, provided it is greater than the numerical precision of the computer processor.

Electric Feld for Three ECT Electrode Configurations
The MRI scan of a patient (39-year-old male) with bipolar disorder was acquired at Neuroscience Research Australia following an ECT session. The patient provided an informed consent for study participation, which also received ethics approval by the University of New South Wales. A T1-weighted 3 T head scan was obtained along with diffusion tensor imaging in 32 directions. The head scan was truncated at the chin, and the voxel size was 1 mm in all directions. The head was segmented into multiple domains corresponding to the tissues listed in Table 2 Three common ECT electrode configurations were simulated: bifrontal (BF), bitemporal (BT) and right unilateral (RUL) as depicted in Fig. 2.7. An electrical ECT stimulus was applied using an isolated circular boundary defined to be 5 cm diameter. This boundary was supplied with a current density J at the anode and -J at the cathode. The stimulus current was set at 800 mA in all electrode configurations. The lower neck boundary was set to a distributed impedance as in Eq. (2.7), with d s set to 5 cm from the ground reference, whilst σ neck was set to the scalp conductivity.
The analysis was performed on regions of the brain associated with emotional responses, directing attention, memory, verbal and learning skills, among others [37], namely the cingulate gyri, parahippocampal gyri, subcallosal gyri, amygdala, inferior frontal gyri, hippocampus and middle frontal gyri. These ROIs were extracted automatically using BrainSuite's labelling tool and BCI-DNI atlas. This approach allows specific comparison in the particular brain region rather than qualitative observation in the whole brain.
Average  E was calculated for each ROI using Eq. (2.5) and the results for each electrode configuration were compared. The results in Fig. 2.8 indicate that the right middle frontal gyrus displayed the largest average E-field for all three electrode configurations. RUL produced a maximum E-field in all right sides of the analysed locations as well as the left cingulate gyrus. Otherwise, the maximum average E of the other left-brain regions was achieved via the BT configuration. The resulting electric field magnitude,  E , is displayed in Fig. 2.9. The electric field in the BF and BT configurations showed a more symmetric electric field profile, whilst RUL displayed a larger  E over the right brain hemisphere relative to the left. The BF configuration resulted in higher  E at the frontal lobe regions of the brain only. On the other hand, the  E distribution in the BT configuration encompassed both the frontal and temporal lobe regions.
The BT configuration has been noted to be the most efficacious among the three electrode placements using the least amount of stimulus dose; however, it is also associated with a larger rate of cognitive side effects [38]. These effects are possibly due to the larger brain area affected by the electrical field as shown in Figs. 2.8 and 2.9. The RUL configuration is known to result in less verbal impairment [39], possibly due to the lesser impact on the speech area at the left inferior frontal region. Nonetheless, the RUL configuration typically requires a suprathreshold stimulus (exceeding seizure threshold) to be effective [39]. The BF configuration, on the other hand, is also known to result in less cognitive impairment than BT, but nonetheless using lesser stimulus dosage than RUL [40]. As such, electrode placement sites and stimulus dosage can affect the outcome of ECT, which can be investigated further using the modelling framework described in this chapter.
Modelling results such as those shown here can provide additional understanding of clinical ECT findings in a given patient. The electric field, voltage and current data obtained from the simulation can be used to infer ECT impact on the brain. Nevertheless, these results are preliminary and should be treated with caution since results may differ from subject to subject. However, such results can be used to associate brain regions affected by ECT with patient responses following the treatment, providing for the future development of safer and more effective ECT protocols.

Model Extensions
The modelling framework described in this chapter only addresses the computation of average electric field within different brain regions. Other possibilities for data analysis include: focality analysis by masking regions with electric field magnitude below a specified threshold [41], reconstruction of binarised subtraction maps for direct comparison of stimulus effects among different electrode placements [3], and analysis of heating during brain stimulation by incorporating a bioheat equation [42].
Moreover, the framework described here only simulates the head as a passive volume conductor under electrical stimulus, disregarding the complex excitable tissue properties of brain neurons. Nevertheless, the underlying brain activity is still Fig. 2.9 Distribution of electric field magnitude in V/m on the grey matter surface. The maximum electric field is capped to 200 V/m for image clarity. The labels on the first row indicate the image orientation as follows: A anterior, P posterior, R right and L left poorly understood. The neuronal tissue itself acts as an internal source of current, which may disrupt the externally applied stimulus. A review by Ye and Steiger summarizes the evidence for such phenomena from various experimental and simulation studies [43]. Brain neurons are also interconnected to more distal regions due to the existence of neural tracts, which connect different regions of the brain. Thus, excitation of one region will propagate along the axon [37].
Several modelling studies have addressed the issue of brain activation. A recent study by Riel et al. [44] performed a preliminary activation analysis along the white matter fibre tracts using a volume conductor model. The fibre tracts were extracted from a DW-MRI image using constrained spherical deconvolution, and electrical potentials were subsequently interpolated along the fibre tracts for calculating the activation function [44]. Bai et al. [36] presented a finite element (FE) whole-head model incorporating Hodgkin-Huxley-based continuum excitable neural descriptions in the brain, which was able to simulate the dynamic changes of brain activation directly elicited by ECT, allowing investigation of parameters such as pulse duration [36]. Nevertheless, the computation was rather lengthy. In addition, the intracellular potential in the model was assumed to be resistively tied to a remote fixed potential, whose physiological meaning was difficult to interpret. This constraint did not allow for the spread of excitation through neural networks in the brain.
McIntyre et al. have, over the years, introduced a representation of white matter (WM) fibres in the vicinity of the subthalamic nucleus (STN), combined with a volume conductor model of deep brain stimulation (DBS) [45,46]. After the electric potential induced by a DBS device was calculated by the FE solver, the timedependent transmembrane potential was solved in NEURON software using a Hodgkin-Huxley-type model based on the interpolated potential distribution along the length of each axon. The model was able to predict activation in STN neurons and internal capsule fibres, and the degree of activation matched well with animal experimental data [45].

Subject-Specific Tissue Conductivity
The electrical conductivities of head tissues in most modelling studies are based on in vitro measurements. It is anticipated that some variation exists due to experimental conditions and sample preparation, let alone any inter-subject differences. However, it would be highly invasive to perform in vivo measurements of electrical conductivity, especially within brain structures.
Recent work by Fernández-Corazza et al. used electrical impedance tomography (EIT) to noninvasively obtain the conductivity of head tissues [4]. They injected small amounts of electrical current through multiple electrode pair configurations. A Laplace equation, as in Eq. (2.3), was then solved for every electrode pair. Afterwards, an inverse problem was solved to optimise the electrical conductivity using Newton's optimisation method to fit the Laplace equations to EIT experimental data. It was observed that the accuracy of this method also depends on the accuracy of the skull segmentation. Their study showed that an over-smoothed and compact skull geometry, with closed foramen, can overestimate the conductivity by almost 30% relative to a more accurate skull segmentation.

Conclusion
Computational head models and FE simulation provide additional insights into understanding regions of the brain affected by ECT and other transcranial stimulation techniques by using metrics such as the electric field distribution, which are difficult to obtain by direct experimental measurement. Furthermore, these head models can serve as a tool for testing novel ECT protocols prior to animal and clinical studies.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.