Electromagnetic Modeling of Rat’s Head: Comparison of Formulations and Models

Reliable inverse imaging of source currents in rat’s brain requires sufficiently accurate and CPU-time moderate forward models of fields to calibrate inverse solvers. In this paper, we compare different mathematical formulations of the electromagnetic problem related to the analysis of brain waves (static, quasi-static, full-wave) and various meshes differing in the density, the type and the geometrical accuracy. A sufficiently accurate model of brain waves is then completed by the cerebrospinal fluid and the skull. The resultant composite model of rat’s head with properly set electrical parameters has to be calibrated by the outputs of measurements. That way, a realistic electromagnetic model of the head of a live rat can be obtained.


Introduction
In the open literature, methods for solving both the forward problem and the inverse one in electroencephalography (EEG) are described for human beings in detail. Current densities excited by equivalent source dipoles are expressed by Poisson's equation with full-tensor anisotropic conductivity. The human skull is approximated by a three-shell concentric spherical model (the skin, the skull, and the brain). For layered spherical anisotropic volume conductors, semi-analytical and numerical solutions are available.
However, neither forward models nor inverse solvers have been properly calibrated and tested in human beings since in-vivo experiments exploiting deep electrodes are not allowed. Thus, all electromagnetic quantities in human beings have been measured indirectly.
In order to implement direct measurements, experiments can be carried out on rats. Nevertheless, there are significant differences between the brain of a human and a brain of a rat:  Shape and size of brains differ significantly. Due to shape and size limitations, animal brains can be accessed from the top of the skull only.
 Location of source currents in brains is considered differently. In the EEG studies of the human brain, currents are assumed to flow in the cortex mainly. Subcortical sources are ignored because of limited sensitivity of EEG to the signals they generate. In case of the animal brain, the whole volume of the brain can be investigated.
 Electromagnetic parameters of live brain tissues have not been determined with a sufficient accuracy yet. Therefore, numerical models and phantoms have to be calibrated.
 The animal head is covered by muscles. Since muscles distort EEG signals, measurements have to be performed on a shaved skin or on the brain surface directly.
The above-given reasons are the main motivation for the research of brain waves propagating in the head of a rat.
Brain waves can be measured on rat's head surface in form of electrical potentials [1]. Localization of brain wave sources requires a forward model and an inverse solver [2]. The forward model is aimed to calculate the potentials at the surface electrodes from known source parameters.
In the matrix form, the forward model can be expressed as [3] B = L  J (1) where B is the matrix of simultaneously acquired time series of potentials measured in independent channels, J is the current matrix and L is the so-called lead-field matrix.
The lead-field matrix is essential for solving the inverse problem The solution is asked to determine wave sources J in the brain from electrical potentials measured on the surface B.
The lead-field matrix can be evaluated by several simulations of the forward model with variating positions of excitation sources. The paper is focused on the forward model.
At present, most models of brain waves are based on the quasi-static formulation. Thus, the relationship between the brain current sources and corresponding potentials on the brain surface is given by Poisson's equation:  In the review paper [4], basic principles of forward modeling are given. Attention is turned to the generation of EEG, Poisson's equation, a three-shell spherical head model and proper application of numerical methods (BEM, FEM, FDM), especially  In [5], the conductivity profile in the somatosensory barrel cortex of the Wistar rat is evaluated. First, boundaries of a six-layer barrel cortex are determined as spherical ones using fluorescent Nissl staining images. After injecting the current and recording corresponding potentials by electrodes, the analytical formula in a non-linear optimization method is used to estimate the conductivity profile in a spherical model (radii identified with fluorescent Nissl staining images).
In [5], the influence of the estimated conductivity profile to the current source density (CSD) is also demonstrated when forward modeling is included.
 In [6], the extensive description of measuring techniques of the neural activity is provided. The most important part relates to the EEG and MEG forward modeling. Additionally, the intracellular and extracellular recording, CSD and the most popular techniques used to measure neuronal activity (functional magnetic resonance imaging, calcium imaging, etc.) are reviewed.
In majority of today's models, potentials are calculated by the boundary element method (BEM) on the surface only. Hence, no anisotropy of brain tissues can be considered. Moreover, the models cannot include potentials inside the head when stimulating electrodes are used, for example. Therefore the finite integration technique (FIT) was used for simulations presented in the paper.
For a human head, an analytical approximate model consisting of three concentric spheres is known [4]. Since the head of the rat cannot be approximated by a sphere sufficiently [7][8][9][10], more realistic approximations have to be used. Such approximations have not been described in the open literature sufficiently yet.
Furthermore, a recent paper [11] on the time-domain full-wave formulation of brain waves showed significantly different results compared to quasi-static approximations of a sphere and a slab.
Due to these reasons, we decided to compare:  Outputs of all solvers of Maxwell equations available in the CST Studio Suite [12] with an analytical solution of electric field intensity excited by a dipole inside a spherical model of a head. The comparison of results should answer the question what solver should be selected to obtain reliable results with minimum CPU costs.
 Simulations of rat's brain model of a realistic shape [7][8][9][10] and a spherical one. The comparison should show dependencies that cannot be clearly provided by a spherical model (sphere is an insufficient approximation of the shape of rat's brain).
 Rat's brain models of a different geometrical accuracy and a mesh density. The comparison of results should answer the question what geometrical accuracy of the model is needed to calibrate the imaging procedure.
 The isolated brain and the complete head. The brain model was embedded into the cerebrospinal fluid and the skull. The fluid and the skull differ in the electrical conductivity and the relative permittivity. The comparison of results should answer the question what way the fluid and the skull transform the electric field intensity measured on the brain surface.
Preliminary results of the above-specified research were discussed in [2], and outputs of discussions were reflected in the paper.

Formulation of Maxwell Equations
The brain can be understood as a volumetric source of a current produced by populations of neurons. Various frequencies of electromagnetic waves can be observed due to the complex behavior of the brain. Generally, many studies on neural oscillations conclude that deeper parts of the brain excite electromagnetic waves with rather lower frequencies and the frequency rises in the evolutionary younger areas closer to the brain surface [13]. However, appropriate mapping of frequencies is not available due to the complexity and inter-individual variability of the brain.
Brain waves propagate at frequencies from 0.1 Hz to 100 Hz. However, the power of the electroencephalograph (EEG) signal decreases with the increasing frequency by the factor of 1/f [14]. Thus, rat's brain was decided to be simulated at the frequency 10 Hz which corresponds to the known peak in the spectral power of a human EEG. of rat's brain in the CST Studio Suite [12]. The models resulting from these formulations differ in the accuracy and CPU-time demands. The models can be characterized as follows:  The static model represents the highest simplification of Maxwell equations. The static model is related to two solvers of the CST Studio Suite [12], to the electrostatics field solver (field excited by stationary charges) and the stationary current solver (field excited by direct currents). Brain waves can be characterized by the spatial distribution of the electric potential or the electric field intensity.
 The electro-quasi-static model neglects displacement currents, and brain waves are represented by the spatial distribution of the field intensity.
 Full-wave models can be formulated in the frequency domain (a steady state at a single harmonic frequency) and the time domain (a transient state with an arbitrary time response). In both the domains, brain waves can be displayed in form of time-varying spatial distributions of field intensities.
In order to develop the full-wave time-domain model, the frequency scaling method [11] was applied since the lowfrequency breakdown can occur, and the solution of Maxwell equations can diverge.
The frequency scaling method [11] is based on scaling the relative permittivity ε rsc of the brain tissue  r rsc r sc Here, f denotes the original frequency (i.e., the frequency of brain waves) and f sc corresponds to the frequency of simulation. Thanks to this scaling, the simulation can be shifted from Hz to MHz. If brain waves are simulated at higher frequencies, then the low-frequency breakdown is overcome, and the simulation is accelerated.
After the simulation, the resulting quantities have to be recalculated according to the following equation [11]: Here, X represents the observed quantity (intensity of electric or magnetic field) at the original frequency f, and X sc corresponds to the quantity of the simulation at f sc .
After removing duplicated surfaces and edges of elements, correcting unevenness, protrusions, sharp edges and modifying the orientation of mesh elements, a new mesh can be generated to preserve anatomy of rat's brain. Brain models were resized to the length l = 27 mm, the width w = 15 mm and the height h = 10 mm which correspond to approximate dimensions of the brain of an adult rat. Two brain models differing in the mesh density were generated (see Fig. 2). The bold curve in the center of models represents the trajectory followed when electric field intensity E on the brain surface is measured. The coordinate p 1 runs from left to right.
To compare the model of the isolated brain and the complete head, the brain model #1 was embedded into the cerebrospinal fluid and the skull. Since the brain model #1 and the skull model come from different MRI scans of rat's head, the brain was proportionally scaled to fit the cavity in the skull. The free space between the brain and the skull was filled in by the cerebrospinal fluid. After that, the brain model #1, the cerebrospinal fluid and the skull were converted to voxels. Thanks to the conversion, the brain, the fluid and the skull can be aligned (potential gaps can be eliminated), and simulations are accelerated.
In Fig. 3, the composite model of the complete head is shown. The model consists of the brain, the fluid and the skull converted into voxels with the length of the edge  l e = 0.7 mm. The bottom picture in Fig. 3 proves that all the models are properly aligned thanks to the conversion to voxels.
In order to validate simulation results, a simplified spherical model of rat's brain was created. The model consisted of a homogeneous sphere with radius r = 50 mm placed at the origin of coordinates.

Simulation Setup
The group of active neurons can be modelled as a current dipole [4], [6], [19] which excites an electrical field propagating to the surface of rat's brain. The dipole has the total length l = 1 mm with the gap g = 0.1 mm between arms and the diameter of arms d = 0.1 mm.
In simulations with the sphere, the horizontal dipole in origin of coordinates played the role of a source. In all other simulations, fields were excited by the vertical dipole in the center of the brain. Sources excited the dipole moment m = 1 Am. Outputs of all solvers were converted to the distribution of electric field intensity.
Electrical parameters of rat's tissues were approximated by the corresponding parameters of human tissues [4], [20], [21]:  All tissues and formulations considered the relative permeability µ r = 1.
 The electric conductivity was  = 0.33 S/m for the brain and the sphere,  = 1.79 S/m for the cerebrospinal fluid, and  = 0.0174 S/m for the skull.  The relative permittivity was ε r = 1 for the sphere, ε r = 4.07·10 7 for the brain, ε r = 1.09·10 2 for the cerebrospinal fluid, and ε r = 5.51·10 4 for the skull.
The frequency tends to zero in the static formulation. Since relative permittivity is unknown at zero frequency, the values corresponding to 10 Hz were used.
Appling the frequency scaling (3) in the full-wave time-domain formulation, the simulation was transferred from 10 Hz to 10 MHz and the relative permittivity was shifted by six orders.
Due to the complexity and inter-individual variability of the brain, general value of the relative permittivity can be hardly determined. From the electrical viewpoint, the gray matter is the most significant part since neural activities are formed there [4]. In the simulations with the model #1, the model #2 and the complete model, the relative permittivity of the gray matter was therefore used.

Spherical Model
In Fig. 4, the analytical solution of electric field intensity inside the sphere is compared with outputs of solvers available in CST Studio Suite [12].
The electric field intensity E was measured in plane of the dipole along x axis. The results indicate the minimum  difference between the electro-quasi-static solver and the analytical solution. On the other hand, the maximum difference appears in case of the static model. Table 2 compares CPU-time demands and percentage deviations from the analytical solution for all solvers available in the CST Studio Suite. The comparison shows that outputs of the full-wave frequency-domain solver and the electro-quasi-static solver differ from the analytical results similarly, but CPU-time demands of the electro-quasi-static solver are 30× lower approximately.
The full-wave time-domain solver shows a higher error than the frequency-domain one due to the application of frequency scaling. The relative error of the full-wave time-domain solver corresponds to the error of the stationary current solver. Figure 5 shows the influence of the number of tetrahedrons used for meshing. For this comparison, the electroquasi-static solver and the brain model #1 were used. The mesh density is specified by the number of cells per the maximum edge of the model box. The meshing procedure is based on the surface model which is identical for all simulations.

Isolated Brain
The comparison shows that the model with 94 084 tetrahedrons provides results similar to meshes with highest numbers of elements. Therefore, the mesh consisting of 94 084 elements can be considered as sufficiently dense. Table 3 summarizes influence of the number of mesh elements on the accuracy and CPU-time demands of simulations. The simulation with the highest number of tetrahedrons is used as a reference for relative error calculation. The results show that the model with the highest number of tetrahedrons needs the highest CPU time. Since the relative error of the model with 94 084 elements is acceptable, this number of tetrahedrons can be considered to be sufficient for forward modeling. Figure 6 compares simulations by the electro-quasistatic solver for models #1 and #2. The discrepancy at p 1 = 6 mm and p 1 = 26 mm is caused by transitions from olfactory bulbs to the central part of the brain, and from the central part of the brain to cerebellum, respectively (see Fig. 2). These transitions are represented by grooves in 3D   models. The surface of the brain is closer to the dipole, and the electric field intensity is higher there. Since the discrepancy is small compared to the size of grooves, the rough brain model can be considered to be sufficient for forward modeling.
Since rat's brain cannot be approximated by a sphere as the human brain, more realistic models of rat's brain have to be used. Figure 7 shows electric field intensity along the trajectory p 1 computed by all the solvers for the simplified rough model (#1) in the top chart, and for the refined model (#2) in the bottom chart. The measuring trajectory p 1 corresponds to the bold curve on the surface of brain models in Fig. 2 and is measured from left to right.
Comparison of both the dependencies shows the following:  The e-static solver and the stationary current one give comparable results both for the rough model #1 and the fine model #2.
Simulation of the spherical model shows a lower relative error in case of the stationary current solver than the e-static one.
 The electro-quasi-static solver and the full-wave frequency domain one give comparable results for both the models.
 Accuracy of the full-wave time domain solver is influenced by the model geometry. For the fine model #2, the accuracy corresponds to the e-static solver and the stationary current one along the entire length of the measuring trajectory p 1 . For the rough model #1, the accuracy is closer to the electro-quasistatic solver in the middle of the measuring trajectory p 1 ∈ (6 -25) mm. In the rest of the measuring trajectory p 1 , the accuracy corresponds to the e-static solver.
Error of the full-wave time-domain solver is caused by the frequency scaling [11] from 10 Hz to 10 MHz. At this frequency, the wavelength is much smaller and the interaction with the model is stronger.
Since both the model #1 and the model #2 is homogeneous and isotropic, higher values of electric field intensity E correspond to active areas in a shorter distance. Therefore, the highest values of electric field intensity E can be measured in the central part of the measuring trajectory p 1 which is in the shortest distance from the dipole located in the center of the brain. Obviously, simulations with realistic 3D models of rat's brain cannot correspond to simulations with the spherical model.
The electro-quasi-static solver [12] can reach the best accuracy within the lowest CPU time, but does not allow to include lumped elements into simulation which is useful when expanding the forward model (EEG measurement elements, RLC circuits simulating coupling between cables, etc.). Fortunately, inclusion of lumped elements is supported by the full-wave frequency-domain solver which reaches a comparable relative error, and can be an alternative when the CPU time is not the priority.
An overview of CPU-time demands and relative errors of different solvers is listed in Tab. 4. The error is related to outputs of the electro-quasi-static solver.
The electro-quasi static solver requires a similar CPU time compared to the full-wave frequency-domain solver when the fine model #2 is simulated. In case of the rough model #1, CPU-time requirements of the electro-quasi static solver are significantly lower. Thus, the full-wave frequency-domain solver is CPU-time effective when complexity of a 3D model is higher (e.g., when a skull, a cerebrospinal fluid, measuring electrodes, anisotropic properties are added).  The full-wave time-domain solver shows much higher CPU-time demands for the fine model #2 in comparison with the rough model #1. Considering both the accuracy and CPU time demands, the full-wave time-domain solver seems to be inappropriate for complex forward modeling. Figures 8 and 9 show the distribution of computed electric field intensity both for the model #1 and the model #2. On the top, electric field distribution on the surface of the brain is depicted. Other figures visualize field distribution in the coronal cut, the axial cut and the sagittal cut passing the center of the brain. The electro-quasi-static solver interacts with the brain model properly thanks to the long wavelength of the excitation signal and small electric dimensions of the model. Using the frequency scaling method [11] in conjunction with the full-wave time-domain solver, the frequency is shifted from 10 Hz to 10 MHz, the wavelength is shortened and electrical dimensions of the model are increased. Thus, the fullwave time-domain solver is not recommended to be used. Fig. 9. Electric field distribution in model #2 evaluated by electro-quasi-static solver. From the top: the brain surface, the coronal cut, the axial cut and the sagittal cut.

Complete Head
The model #1 of rat's brain was converted to voxels, and completed by the cerebrospinal fluid and the skull. The measuring trajectory p 2 (the blue lines in Fig. 3) are measured from left to right (the horizontal axis in Fig. 10 and Fig. 11). Due to the stairway surface of voxel models, the measuring trajectory was created on the top of models to follow the planar surface. That way, continuous dependencies without quick variations caused by sharp edges of voxels were obtained. Moreover, the measuring trajectory on the brain, the fluid and the skull can be of the same length.
In Fig. 10, the influence of the cerebrospinal fluid and the skull is shown. For the simulation of the brain model #1 completed by the fluid and the skull, the electroquasi-static solver was used to compute electric field intensity along the trajectory p 2 . Comparison shows that the field distribution on the pure brain differs both in the magnitude and the shape. The attenuation of electric field intensity E is caused by transitions from the brain to the cerebrospinal fluid and from the fluid to the skull.
Electric field intensity on the pure brain reaches the maximum magnitude in the smallest distance from the excitation source. Increasing the distance from the source, the electric field intensity gradually attenuates. Nevertheless, the addition of the cerebrospinal fluid and the skull causes irregular distribution of electric field intensity E along the measuring trajectory p 2 .
The results show two parts of measuring trajectory p 2 where the electric field intensity E is higher for the complete model of the head. However, the results for the rest of the measuring trajectory p 2 indicate smaller values of electric field intensity E. Moreover, the electric field intensity E is attenuated faster when the distance increases from the maximum value compared with pure brain. Thus, the inclusion of more parts of rat's head enhances the accuracy of the forward modeling.
The complete head was simulated by the electroquasi-static solver. The measuring trajectory p 2 was on the surface of separate models (see blue lines in Fig. 3). The position was the same for x, z coordinates and differed in y coordinate. A different y position corresponds to separate models (the brain, the cerebrospinal fluid, and the skull).
In Fig. 11, the distribution of electric field intensity E on surfaces of separate models is irregular due to the interface among models. Hence, the measuring position in rat's head is important for forward modeling and EEG measure-  ment. In real EEG measurements, a small shift in the depth of measuring electrodes can give different results. Figure 12 shows electric field distribution evaluated by the electro-quasi-static solver for the complete head on the surface of the brain, in the coronal cut, the axial cut and the sagittal cut passing the center of the head. In the picture of the field distribution on the brain (the top picture), the cerebrospinal fluid and the skull are hidden. The result indicates influence of separate head models to the distribution of electric field intensity.

Conclusions
In the paper, electromagnetic models of rat's brain and rat's head were described. Attention was turned to the discretization of the brain, the selection of the appropriate solver, and the composition of the whole head of the brain, the cerebrospinal fluid and the skull.
When applied to the spherical model of rat's head, comparison of CST solvers shows that the electro-quasistatic solver produces the minimum relative error compared to the analytical solution. On the other hand, the e-static solver shows the highest inaccuracy compared to the analytical solution. The frequency scaling method applied in the full-wave time domain method increases the relative error and produces results comparable to the stationary current solver. The full-wave frequency domain solver achieves the relative error comparable to the electro-quasistatic solver, but CPU-time demands are higher.
When changing the number of mesh elements, almost the same dependencies of electrical fields are obtained. The model defined by 94 084 tetrahedrons results in the same dependency as the model with the highest number of tetrahedrons. Nevertheless, the CPU-time demands are much higher for the model with the highest number of tetrahedrons. Thus, the smaller number of tetrahedrons is sufficient for forward modeling.
Since rat's brain cannot be approximated by a sphere, more realistic models of rat's brain were developed. The same electric-field dependencies were obtained when the e-static solver and the stationary current one were used. Nevertheless, the stationary current solver showed much smaller relative error compared to the e-static solver when the spherical brain was simulated. The stationary current solver is not suitable for realistic forward modeling.
Accuracy of the full-wave time-domain solver is related to the accuracy of the model geometry:  In case of the fine model #2, the accuracy is comparable with the e-static solver and the stationary current one.
 In case of the rough model #1, the accuracy approaches the electro-quasi-static solver.
The error of the full-wave time-domain solver is caused by the frequency scaling method [11] which shifts the frequency of analysis from 10 Hz to 10 MHz. The higher frequency corresponds to the smaller wavelength, and the interaction of waves with the model is different.
The best results from the viewpoint of accuracy and CPU-time demands were reached by the electro-quasistatic solver which does not support simulations of lumped elements. Fortunately, simulations of lumped elements are supported by the full-wave frequency-domain solver with a similar relative error and higher CPU-time demands.
Comparison of all the solvers from the viewpoint of CPU-time demands showed that the e-static solver and the stationary current solver can save the CPU time compared to other solvers in case of the rough model #1. For the fine model #2, the CPU time demands were similar for all the solvers except of the full-wave time-domain one.
Comparison of all solvers shows that:  The e-static solver and the stationary current one excel in CPU-time demands but fail in accuracy.
 The full-wave time-domain solver is not suitable for forward modeling if both the accuracy and CPU-time demands are considered at the same time.
 Both the electro-quasi-static solver and the full-wave frequency-domain solver excel in accuracy. When simulating simpler models, the full-wave frequencydomain solver requests more CPU time.
In order to create the model of the whole head, the rough brain model converted into voxels (the model #1) was completed by the cerebrospinal fluid and the skull. Thanks to the voxel conversion, the perfect alignment of the brain, the fluid and the scull was ensured, and the simulation was accelerated.
Comparison of the separate brain (converted to voxels) and the complete head indicates that dependencies for the separate brain differ in the magnitude and the shape. If the cerebrospinal fluid and the skull are added, the attenuation of the calculated signal can be observed. Moreover, the irregular distribution of electric field intensity was shown. If more parts of rat's head are included, a better accuracy of forward modeling is reached.
For the complete head, the irregular dependencies of the electric field intensity were obtained for the brain, the fluid and the skull surfaces. Therefore, the small change in depth can give a different value of electric field intensity. This irregular distribution is caused by interfaces of models with different permittivity and electric conductivity. Spatial distributions of electric field intensity for the model #1, the model #2 and the complete model were visualized on the surfaces, and in the cut passing the center of the head (the coronal plane, the axial plane and the sagittal plane were considered). The result indicates the influence of separate parts of the head to the distribution of electric field intensity.
In near future, the numerical model of rat's head will be validated experimentally. Electrical potentials will be measured on the brain surface of a live rat in 13 different points. The numerical model of the brain will be calibrated to reach a reasonable agreement. The process of the calibration should indicate whether an isotropic model of the brain is sufficient for forward modeling, or anisotropy of the brain tissue is necessary to be considered.