Computational Models of Brain Stimulation with Tractography Analysis

Computational human head models have been used in studies of brain stimulation. These models have been able to provide useful information that can’t be acquired or difficult to acquire from experimental or imaging studies. However, most of these models are purely volume conductor models that overlooked the electric excitability of axons in the white matter of the brain. We hereby combined a finite element (FE) model of electroconvulsive therapy (ECT) with a whole-brain tractography analysis as well as the cable theory of neuronal excitation. We have reconstructed a whole-brain tractogram with 2000 neural fibres from diffusion-weighted magnetic resonance scans and extracted the information on electrical potential from the FE ECT model of the same head. Two different electrode placements and three different white matter conductivity settings were simulated and compared. We calculated the electric field and second spatial derivatives of the electrical potential along the fibre direction, which describes the activating function for homogenous axons, and investigated sensitive regions of white matter activation. Models with anisotropic white matter conductivity yielded the most distinctive electric field and activating function distribution. Activation was most likely to appear in regions between the electrodes where the electric potential gradient is most pronounced.

models for a better understanding of the mechanisms underlying brain stimulation, a passive model unable to mimic active membrane dynamics is of limited utility. Bai et al. presented a finite element (FE) whole head model with the incorporation of a Hodgkin-Huxley-based continuum excitable neural description in the brain, which was able to simulate the dynamic changes of brain activation directly elicited by ECT [2]. 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 find in reality. This constraint also led to the missing of the spread of excitation through neural networks in the brain (seizure). McIntyre's group has over the years introduced a representation of the white matter (WM) fibres in the vicinity of the subthalamic nucleus (STN), combined with a volume conductor model of deep brain stimulation (DBS) [3,4]. After the electric potential induced by the DBS device was calculated by the FE solver, the time-dependent transmembrane potential was solved in the neuron using a Hodgkin-Huxley-type model based on the interpolated potential distribution along the length of each axon. The model was able to predict the activation in STN neurons and internal capsule fibres, and the degree of activation matched well with animal experimental data [3]. This chapter describes the steps necessary to implement computational modelling of the human head, including a white matter tractography analysis to determine the voltage, electric field and activation function distribution for two prescribed setups of ECT. Neural activation is analysed by the neural activation function, which is in turn approximated by the second spatial derivative of the electrical potential, assuming that the axonal diameter, axoplasmic resistivity and capacity were constant.

Methods
To evaluate the electric field and the neural activation function along neural fibres, an FE volume conductor model of electrical brain stimulation (ECT in this case) was combined with a white matter fibre tractography model of the brain. Figure 1 details all the necessary steps for a combined tractography analysis of electrical brain stimulation.

Image Preprocessing
The FE volume conductor model and the white matter tractography model were both reconstructed from magnetic resonance image (MRI) data from the same subject. MRI data were taken from subject MGH1010 (gender: female, age: 25-29) provided by the NIH Human Connectome Project (HCP) [5]. The dataset comprised a structural T1-weighted scan, a high-resolution T2-weighted scan and a set of highresolution, high b-value diffusion-weighted (DWI) scans. These MR scans provided by HCP had been preprocessed, and therefore image distortions/artefacts currents caused by eddy, gradient nonlinearities and motion had been removed [6]. To ensure the same frame of reference for the individually acquired image datasets, it is necessary to perform image registration. The structural image data with a resolution of 1 mm was resampled to a resolution of 1.5 mm using trilinear interpolation in order to match the dimension of DWI data. Subsequently, a rigid-body registration with six degrees of freedom was performed between the fixed DWI data and the floating T1-weighted image. Since the transformation was performed on the structural data (aligned to the DWI data), no further processing of auxiliary diffusion data (gradient tables of DWI scans) was necessary. These operations were implemented using FLIRT (FMRIB's Linear Image Registration Tool) provided in the FMRIB Software Library (FSL) which is developed by the Oxford Centre of Functional MRI of the Brain (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT) [7-9].

White Matter Fibre Tractography
White matter consists of myelinated neuronal tracts, interconnecting different parts of the brain. This neural fibre network causes a highly anisotropic diffusion behaviour of water molecules inside the brain. Diffusion is higher along the fibres than in the transverse direction [10]. This characteristic is utilised by DWI scanning, which Flowchart describing the workflow to generate a finite element volume conductor model and a white matter fibre tractography of the head. Information on the electrical field and the activation function can be computed by numerical approximation of the spatial derivatives of the voltage along the fibre, which in turn can be extracted by combining both finite element head model and the tractogram applies a magnetic gradient in different directions during image acquisition and thus encodes the sensitivity of diffusion to each gradient direction in the image data. Similar anisotropic behaviour is found for the electrical conductivity of white matter fibres, and a linear relationship between electrical conductivity and water diffusion has been experimentally validated [11,12]. Therefore using anatomically constrained tractography, we extracted a geometrical model of neural fibres from DWI data, which provided information about the trajectories of white matter fibres. Structural T1-weighted MRI data containing information on the structural composition of the brain was used to determine the start and end of each fibre.

Image Segmentation
In order to adopt anatomical constraints for tractography, the structural data was firstly segmented. T1-weighted image data is typically used as it provides a good contrast between head structures, especially between grey and white matter [13]. Image segmentation was performed in two steps: brain extraction and brain tissue segmentation. The brain extraction procedure was implemented in the FSL Brain Extraction Tool (BET) (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/BET/UserGuide). An intensity histogram of the structural images was used to find the threshold to differentiate the brain-and non-brain regions. Triangular tessellation of a spherical surface was initialised and iteratively deformedbased on the brain/non-brain intensity thresholdto wrap the whole-brain volume [14]. The extracted brain volume was then segmented into cortical grey matter, subcortical grey matter, white matter, cerebrospinal fluid and pathological tissue. The procedure was executed in the FLS automated segmentation tools: FAST (https://fsl.fmrib.ox.ac.uk/fsl/ fslwiki/FAST) and FIRST (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FIRST). FSL FAST makes use of a hidden Markov random field method to categorise grey matter, white matter and CSF [15]. Due to poor contrast, an intensity-based segmentation procedure failed to segment subcortical grey matter, but FSL FIRST, which utilises a priori knowledge on the subcortical grey matter, improved segmentation results [16].

Fibre Orientation Distribution
The basis for every white matter tractography is the identification of fibre orientations from DWI data. The fibre orientation distribution (FOD) holds information on the probability of certain fibre directions within each voxel. The measured signal S (θ, ϕ) at a direction (θ, ϕ) is linked to FOD via the axially symmetric response function R(θ) that satisfies where R(θ) models the signal expected for a voxel containing a single, coherently oriented fibre population. FOD was obtained through constrained spherical deconvolution using spherical and rotational harmonics that formed a complete orthonormal basis set of functions over the sphere and the space of pure rotations. A constraint on positive fibre directions was then imposed [17][18][19]. The procedure was performed in the MRtrix3 software package (https://www.mrtrix.org/), which provides tools for the analysis of DWI data [20].

Anatomically Constrained Tractography
FOD estimates the probability that white matter fibres are aligned in a given direction. When deriving a probability density function (PDF) from the FOD, the trajectory of a neural fibre can then be constructed.
Tractography started with seeds being randomly scattered within a defined seeding area. From each seed point, the direction of a step was calculated by randomly sampling the PDF derived from a second-order integration over FODs [21]. The step size was predefined. Instead of stepping along straight-line segments defined by the step size, the connection between two points was realised by an arc of a circle of fixed length (step size), tangent to the current direction of tracking at the current point. The probability of each path was calculated as a joint probability of each infinitesimal step, thus making up the path. This joint probability was approximated by computing the product of the amplitude of the FOD, evaluated at four points (on the arc segment) along the tangent to the path. The sampling of the PDF was achieved using rejection sampling, a simple Monte Carlo method. Streamlines were iteratively generated by stepping along a path determined by the step size, the orientation derived from the PDF and a constraint on the maximum curvature of the streamline. The procedure was repeated until a stopping criterion was met.
Since reconstructed streamlines represent white matter fibres connecting spatially distant areas of grey matter, the start and end points should occur within the grey matter. Some fibres may project to the spinal cord, but fibres should never end within the white matter or fluid-filled regions. This a priori information was utilised to define the seeding area and to derive a termination criterion. Only fibres starting and terminating within cortical or sub-cortical grey matter were accepted. Fibres that ended in fluid-filled regions or inside the white matter were rejected [22].
As the density of reconstructed streamlines does not necessarily represent the true density of anatomical connections, the raw tractogram requires a filter. To find a subset of streamlines that best matched the diffusion signal, a sphericaldeconvolution informed filtering procedure was used [23]. The FOD was proportional to the DWI signal and to the volume of fibres within a volume pixel. The FOD was firstly sampled using a large set of basis direction of a unit hemisphere. The sampled FOD was then segmented into individual lobes representing FOD peak orientations. The tractogram was subsequently mapped to the voxel grid of the DWI scans, and for every voxel, a FOD lobe look-up table was generated. Streamlines within every voxel were assigned to a lobe. For every voxel, the amount of streamlines in a specific orientation (of unit length) needed to match the amplitude of the FOD peak orientations. Therefore, a cost function, based on the FOD lobe amplitudes over all voxels and the total track density, was introduced and minimised by removing individual streamlines from the tractogram. Key characteristics of the tractogram and input settings are listed in Table 1.

Post-Processing
Two more processing steps are necessary before the white matter geometry can be combined with the volume conductor FE model: model reduction and fibre smoothing. The filtered tractogram still contained over 200,000 individual white matter fibres that ranged from 3 to 150 mm. In order to optimise the computation time for a preliminary analysis, the total number of fibres in the tractogram was reduced. Figure 2 shows a white matter tractogram for 1000, 2000, 5000 and 10,000 fibres ranging from 50 to 150 mm. While 1000 streamlines appeared too sparse for a proper visual inspection of the activating function, 2000 and more fibres covered the wholebrain volume in an adequate density. To limit the computational expense, a white matter model of 2000 fibres was selected. The extracted streamlines did not follow a smooth trajectory since the tractography algorithm allowed an abrupt change in direction after every step (within the fibre curvature constraint). To reduce the risk of potential errors and noise due to spatial discontinuity, the trajectory of every streamline was smoothed and subsequently oversampled. A cubic spline was used to fit the streamline data for all three dimensions separately passing through all provided points on the streamline. Every streamline was oversampled by a factor of 10, increasing the number of points and reducing the step size from 0.75 mm to 0.075 mm. This was done using the univariate spline interpolation function provided by the Python SciPy package (https://www.scipy.org/).

Finite Element Analysis of ECT Brain Stimulation
To simulate the potential field introduced to the head by ECT, a volumetric FE model of the head was generated from structural MRI data. This included image segmentation, FE mesh generation, electrical property settings and finally running simulations in a numerical solver. Two different ECT electrode placements were simulated in this study, and three different electric conductivity settings were chosen for the white matter compartment. The FE analysis of ECT brain stimulation was based on Bai et al. [25].

Finite Element Model Reconstruction
The structural T1-weighted MRI data was used to reconstruct an FE head model. After contrast and edge enhancements of the structural data, segmentation of the scalp, skull, paranasal sinuses, cerebrospinal fluid, grey matter and white matter were performed using 3D Slicer, an open-source software platform for medical image processing (https://www.slicer.org/). Each tissue compartment was assigned with a label map using a range derived from a histogram, representing the grey levels of the desired tissue type. The selection of desired tissue was performed by setting the threshold on a single slice, and then a paintbrush was used to manually modify and correct the selection. This procedure was repeated at every second or third slice of the image volume, and a full segmentation was created by using a trilinear interpolation to automatically connect the sparse set of contours. A paintbrush was then used again to further modify the tissue map until the desired accuracy was achieved. As the face of the subject was removed from the images to keep the subject identity anonymous, a smooth surface was used to replace the face of the head model. For every tissue compartment, a surface triangular mesh was generated in 3D Slicer. To increase the mesh quality and smoothness, the surface meshes of the head model were imported into Geomagic Wrap (3D Systems, USA). Non-manifold edges were removed, self-intersecting triangles were split, the edge crease was reduced, spikes were smoothed, and existing holes were repaired. The processed surface meshes were then transferred to ANSYS ICEM CFD (ANSYS, USA) to create a volumetric mesh. Edges at the intersections between compartments and at the desired electrode contact locations were defined, and a tetrahedral volumetric mesh with appropriate meshing and coarsening parameters was generated. To calculate the voltage data, the final volumetric mesh, with approximately half a million tetrahedral elements, was exported to COMSOL Multiphysics (COMSOL Inc., Sweden), a cross-platform FE solver.

Tissue Conductivities
To calculate the potential field induced by ECT, electrical properties need to be assigned to the different tissue compartments of the volumetric head model. All parts of the head model, apart from the white matter compartment, were considered electrically homogeneous and isotropic. Conductivity values assigned to the scalp, the skull, the paranasal sinuses, the CSF and grey matter are listed in Table 2. Three different white matter conductivity settings, labelled with G, GW and anisotropic, were simulated in order to compare its influence on the potential distribution and the derived electric field and activating function: • G: isotropic grey matter conductivity (0.31 S/m) • GW: isotropic conductivity (0.14 S/m) • Anisotropic: anisotropic conductivity based on conductivity tensor [26,27]

White Matter Conductivity Anisotropy
The linear relationship between the electrical conductivity tensor and the water diffusion tensor implies that both share the same eigenvectors [11,12]. The water diffusion tensor was extracted using the diffusion tensor fitting algorithm for DWI data implemented in FSL DTIFIT (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT). The conductivity tensor σ C of a white matter finite element was modelled as where E N is the orthogonal matrix of normalised eigenvectors of the diffusion tensor. σ long C and σ trans C are the eigenvalues parallel and perpendicular to the fibre directions with σ trans C σ long C [27]. The eigenvalues of the conductivity tensor were determined using a volume constraint proposed by Wolters et al. [26]. To ensure that only white matter voxels with a high level of anisotropy were considered for the calculation of the conductivity tensor, the fractional anisotropy (FA) derived from the diffusion tensor was used. Regions with a low FA values (typically FA < 0.45), which suggested a low local anisotropy, were removed from the tensor analysis.

ECT Brain Stimulation Settings
In the field of bioelectromagnetism, since the frequency of bioelectric events is low, biological tissues are typically considered as volume conductors, in which the capacitive and inductive components of the electrical impedance is neglected [28]. Thus, all head compartments were formulated as passive volume conductors, and the electric potential resulting from ECT was obtained by solving Laplace's equation where V is the electric potential and σ C is the electrical conductivity tensor. ECT electrodes were defined as circular boundaries with a radius of 2.5 cm on the scalp. A stimulus of 800 mA was delivered through the stimulating electrode. Two bipolar ECT electrode placements were simulated as Furthermore, a constant current density J n for the stimulating electrode was defined as where S E is the area of the electrode and I S is the applied stimulus current at a constant DC level of 800 mA. All external boundaries were assigned as electric insulators, and a continuous current density across all interior boundaries was assumed.

Model Combination
Having both models set up, the FE volume conductor model, solved for the potential distribution within the head for a given electrode placement, can be combined with the white matter tractography model. For every white matter fibre conductivity setting and ECT electrode position, voltage data along the fibres in the tractogram was exported from the COMSOL models. For fibre coordinates not coinciding with a FE mesh point, trilinear interpolation was used. To investigate neural activation, the neural activating function was evaluated to link the activation to the second spatial derivative of the potential along the fibre [29]. In Eq. (5), r ax is the axonal resistance determined by the fibre diameter and the axoplasmic resistivity, and c m is the capacity of the fibre. The activating function was in this chapter approximated by the second spatial derivative of the electrical potential assuming that r ax and c m were constant. The electric field (EF) E = 2 ∇ V is the negative first spatial derivative of the potential along the fibre. The first spatial derivative ∂V ∂x was approximated by central difference: where k represents the kth node on an individual fibre, V 0 k is the first spatial derivative of V at the kth node and |r k + 1, k | and jr k, k À 1 j are the distances between the (k + 1)th and kth nodes as well as between the kth and (k À 1)th nodes, respectively. The second spatial derivative along the fibre ∂ 2 V ∂x 2 was approximated by the central difference method on the electric field such that 3 Results

White Matter Fibre Tractography Model
The final neural fibre modelas shown in Fig. 4 consisted of 2000 individual fibres with the length between 50 and 150 mm. The spatial sampling of each fibre was 0.075 mm and every fibre respected the anatomical boundary condition of starting and ending within the grey matter. Fibres were manually categorised to commissural fibres (red, connecting both hemispheres), projection fibres (blue, connecting the brain to the spinal cord) and association fibres (green, interconnecting brain areas within the same hemisphere). The fibre of the pons, connecting both hemispheres of the cerebellum (red), marked only a small fraction of fibres. The model contained 58% association fibres, 24% commissural fibres (including fibres of the pons) and 18% projection fibres.

Electric Field and Activating Function for Three White Matter Conductivity Settings
The potential distribution as well as the approximation of the electric field and the activating function along the fibres for a bitemporal and a right unilateral electrode placement is presented in Fig. 5. The comparison of three different conductivity settings for white matter shows that the incorporation of white matter anisotropy resulted in a wider range and a greater non-uniformity of the potential distribution. The non-uniform potential distribution also led to higher amplitudes and a more distinctive distribution of the first and second spatial derivatives when compared to the two isotropic conductivity settings. In addition, the influence on deep brain structures was most pronounced in the models with anisotropic white matter conductivity. Figure 6 shows the electric potential as well as the first and second spatial derivatives of the potential along the fibre for a single commissural fibre in the anisotropic bitemporal ECT model. High-potential values appeared on the left side of the fibre, i.e. the left hemisphere. The fibre started with a steep descent from the right hemisphere towards the corpus callosum where it ran nearly parallel before it took a turn up and down. The symmetric electrode placement resulted in a strong potential gradient from the right hemisphere to the left. Given the trajectory of the fibre, two peaks appeared in the first spatial derivative resulting in a high second spatial derivative at the beginning and the end of the fibre. While crossing the corpus Electric potential as well as the first and second spatial derivatives of the potential along the fibres for bitemporal (left) and right unilateral (right) electrode placements using three different white matter conductivity settings. In 'G' white matter conductivity is assigned the same value as the grey matter (0.31 S/m); in 'GW', the white matter conductivity is assigned an isotropic value of 0.14 S/m; in 'anisotropic', the white matter fibre conductivity is an anisotropic conductivity tensor derived from the fibre orientation distribution callosum, the electric field was nearly constant and thus causing the approximation of the activating function to be close to zero. Figure 7 shows four views of the second spatial derivative for both ECT models with anisotropic white matter conductivity. For a bitemporal electrode placement, association fibres running through the frontal, temporal and parietal lobes of both hemispheres exhibited a higher probability of activation. For right unilateral electrode placement, activation was more likely to initiate from the right hemisphere, among the splenium of corpus callosum and the association fibres that travelled between the posterior part of the frontal lobe and the parietal and temporal lobes.

Discussion
This chapter presented a study that utilised, for the first time, a white matter tractography analysis in an FE model of ECT, which can be applied to models of other types of brain stimulation. In reality, the brain is made up of a tight package of neuronal somas in the grey matter and their extensionsaxonsin the white matter connecting to other neurons. According to the theory of external stimulation of neurons [29], neuronal excitation of a homogenous axon can be indicated by the second spatial derivative of the potential V, i.e. the negative derivative of electric field along the axon. The assumption of a homogenous axon is only valid in white matter. In the grey matter, where cell somas are located, changes in the diameter of the axons to the soma violate this condition, and the precise electroanatomy of the neuron, together with the variation of the extracellular potential V, has to be evaluated. We have therefore combined a tractography analysis with a volume conductor model and derived brain locations in the white matter, where axons are most likely activated by ECT. It should nevertheless be noted that a more accurate analysis of the activation pattern can only be performed by the inclusion of neuronal cable models. Moreover, the addition of cable models can also be used to simulate the propagation of activation, such as seizure. However, as neurons in varied functional areas of the brain have been found to possess different geometrical and ionic channel properties, the incorporation of cable models across the whole brain may be a challenging and computational expensive task. To avoid large numerical errors when the second derivative of V along the path of the axon was calculated, it was necessary to smooth the path of the axon and interpolate potential data carefully onto the fibre trajectory. The purpose of the tractography was to study the distribution of electric signals along the fibres on a large scale. Therefore, it is important to have a correct representation of all types of white matter fibres within the brain. For this preliminary study, however, after a downsampling process to select long fibres, the full variety of neural pathways was not well presented, particularly the short projection fibres. In spite of this, a close inspection of the filtered tractogram revealed a good representation of major fibres connecting the two hemispheres, as well as the cerebrum to the brainstem and the cerebellum. Furthermore, it was confirmed that the fibres did not cross liquid-filled areas and that each fibre started and ended within the grey matter. 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.