Mathematical Modeling of Brain Swelling in Electroencephalography and Magnetoencephalography

: In the present paper, the forward problem of EEG and MEG is discussed, where the head is modeled by a spherical two-shell piecewise-homogeneous conductor with a neuronal current source positioned in the exterior shell area representing the brain tissue, while the interior shell portrays a cerebral edema. We consider constant conductivity, which assumes different values in each compartment, where the expansions of the electric potential and the magnetic ﬁeld are represented via spherical harmonics. Furthermore, we demonstrate the reduction of our analytical results to the single-compartment model while it is shown that the magnetic ﬁeld in the exterior of the conductor is a function only of the dipole moment and its position. Consequently, it does not depend on the inhomogeneity dictated by the interior shell, a fact that veriﬁes the efﬁciency of the model.


Introduction
Human brain activity is known to produce electric and magnetic fields, both in the interior and in the exterior of the head [1][2][3][4][5]. The calculation of these fields considering a given electrochemical source constitutes correspondingly the forward problem of Electroencephalography (EEG) and Magnetoencephalography (MEG), which are two of the most efficient techniques for studying the functional brain. Reference [6] completely answers the fundamental mathematical question of the existence and uniqueness of the representations obtained using these techniques and also covers many other concrete results for special geometric models of the brain, presenting the research of the authors and their groups in the last two decades. The electrochemical source is commonly represented by a point current dipole lying in the interior of the conductive brain tissue, where expressions of the produced magnetic fields in four basic volume conductor shapes can be found [2]. A large amount of noteworthy research has been conducted considering a homogeneous spherical model for the human head [7][8][9] or geometrically much more complicated models [2,10,11], such as the homogeneous ellipsoidal one [12][13][14][15]. Indicatively, in [9], the basic mathematical and physical concepts of the biomagnetic inverse problem are reviewed, wherein the incorporated forward problem has been introduced for both homogeneous and inhomogeneous media. Then, more sophisticated models followed, such as that in [10], in which the computational and practical aspects of a realistically shaped multilayer model for the conductivity geometry of the human head are presented. On the other hand, the complexity increases when the head is geometrically represented by the most general ellipsoidal coordinate system that reflects the complete anisotropy of the three-dimensional space. Toward this direction, the magnetic induction field in the exterior of an ellipsoidally inhomogeneous, four-conducting-layer model of the human head is obtained analytically up to its quadrupole approximation [12], while in [13], the octapolic contribution of the dipolar current to the expansion of the magnetic induction field is provided. Additionally, other models have considered the head as a non-homogenous conductor, in the sense that it is comprised of multiple layers with different electric conductivity [16], representing the cerebrum, the fluid layer, the scull, and the scalp [17][18][19], while others have considered perturbations in specific layers representing tumors or injuries [20,21], in the aim of continuing to advance the understanding of how sensitive the solution of the forward EEG problem is in regard to the geometry of the head.
Brain swelling [22,23] is the accumulation of water in various spaces of the brain (edema), commonly seen by pathologists and clinicians as a common and often nonspecific finding in a wide variety of cerebral disorders in association with tumors, trauma, and infections, as well as with toxic, anoxic, and metabolic disorders. It typically causes compression of the brain tissue and blood vessels within the skull, as well as impaired nerve function. Its diagnosis is commonly confirmed with CT scans and MRI [24,25]. To this end, and in relation to the forward EEG/MEG problem, we propose a piecewisehomogeneous spherical model for the head, which takes into consideration the existence of fluid at the core of the brain. In particular, we assume that the head is a spherical conductor, which consists of an interior conductive sphere representing the accumulated fluid and an exterior concentric spherical cerebrum shell, which contains a current dipole source. A clockwise spherical coordinate system is set appropriately in such a way so as the origin coincides with the center of the physical system. Then, we calculate in closed form the produced electric potential and the magnetic field in the interior and the exterior of the head as functions of the dipole moment and location as well as the fluid core diameter. The complete analytical formulae are given as series expansions in terms of spherical harmonic eigenfunctions, wherein the effect of the source field is incorporated into the expressions. This paper is organized as follows. In Section 2, we present the geometric model and the mathematical formulation of the forward EEG and MEG problems. In Section 3, we evaluate the electric potential in the interior and the exterior of the head, while Section 4 is devoted to the reduction of the results to the single-compartment model, followed by a numerical implementation in Section 5. In Section 6, we provide analytical expressions for the produced magnetic field, while in Section 7, we calculate this magnetic field in the exterior non-conductive space. Finally, we discuss our results and conclude in Section 8.

Statement of the Problem
In the present model, depicted in Figure 1, we consider the human head as a spherical conductor surrounded by the non-conductive space Ω e (air) and formed by the compartments Ω f and Ω c , that are defined by the concentric spherical surfaces S f and S c , corresponding to the radii f and c, respectively. The corresponding electric conductivities σ f and σ c are considered constants with σ f = σ c , whereas the magnetic permeability is considered equal to µ 0 everywhere. The localized electric activity of the cerebrum tissue is represented by an equivalent current dipole located at the fixed point 0 τ inside the region c Ω , which is characterized by the dipole moment Q and the primary current density The localized electric activity of the cerebrum tissue is represented by an equivalent current dipole located at the fixed point τ 0 inside the region Ω c , which is characterized by the dipole moment Q and the primary current density where δ denotes the Dirac measure. This primary current produces an irrotational electric field E, which is represented by the electric potential u, i.e., The electric potential, denoted by u f , u c , and u e in Ω f , Ω c , and Ω e , respectively, satisfies the Poisson equation in the region that includes the source (Ω c ) and the Laplace equation in Ω f and Ω e , i.e., while on the boundary surfaces S f and S c , the electric potential and the normal component of the electric field are continuous In order to obtain a unique solution for the exterior potential problem, we also have the requirement that whereas lim |r|→0 u f (r) = ±∞. (11) Note that, for every τ ∈ Ω f ∪ Ω c ∪ Ω e , the electric potential u can be expressed as where is the Heaviside function. Furthermore, the primary current J p induces the volume currents and in the conductive regions Ω f and Ω c , respectively, resulting in the total current density The total current J t generates the electric field E and the magnetic induction B in both the interior and the exterior space, which satisfy the quasi-static approximation of Maxwell's equations [26] ∇

The Electric Potential
The electric potential solves the system of mixed boundary value problems (3)-(11). This can be split into three separate problems (3), (4), and (5), subject to the boundary conditions (6)-(9) and the limiting conditions (10) and (11), where we have to deal with the solution of Laplace and Poisson equations in spherical coordinates. To this end, we will use the general expansion of the harmonic function u (∆u = 0) where A m n and B m n , for every n ≥ 0 and m = −n, −n + 1, . . . , n − 1, n, are arbitrary constants and Y m n are the complex spherical harmonics which satisfy the orthonormality relation where n ≥ 0, m = −n, −n + 1, . . . , n − 1, n and m n = while P m n denotes the associated Legendre polynomials and Y −m In the first place, Equation (3) represents the harmonic electric potential in the interior region Ω f and, since it refers to an interior BVP that includes the origin, expansion (19) via the condition (11) yields Next, Equation (5) represents the harmonic electric potential in the exterior region Ω e ; hence, as it refers to an exterior problem, utilizing (19) together with (10) gives where the term for n = 0 falls as 1/τ, which is consistent with condition (10).
Nevertheless, the solution of the Poisson equation (4) is not straightforward. Namely, its solution is a superposition of a harmonic function u c,h and a particular solution u c,p that satisfies (4), i.e., where and Equation (26) is a classical Laplace equation that refers to the intermediate region Ω c ; therefore, it has the general solution (19), which is On the other hand, in view of Equation (27), applying the operator Q · ∇ to the two parts of the equation we are led to the particular solution or Then, we utilize the formula where from which we deduce that the fundamental solution of the Laplacian has the following expansion in terms of spherical harmonics where g m for every n ≥ 0, |m|≤ n , while ν depends upon the position τ = (τ, ϑ, ϕ) over the point τ 0 = (τ 0 , ϑ 0 , ϕ 0 ) according to (33). As a result, and bearing in mind that, in general, the operators ∆ and ∇ act on the position vector τ (unless otherwise specified) and considering the identity the particular solution (30), in view of (35), (36), and (37), is written where for every n ≥ 0, |m|≤ n . Therefore, the sought-out solution (25) of the Poisson equation (4), in view of (28) and (38), is written as Subsequently, we proceed to the calculation of the unknown coefficients A m n, f , B m n,e and A m n,c , B m n,c of the expansions (23), (24), and (40), respectively, for every n ≥ 0 and |m|≤ n , utilizing the continuity conditions (6)-(9) and the orthogonality relation (21). In particular, substituting (23) and (40) in (6), for τ = f and in light of (33), we get which, by virtue of (21), leads to the equation for every n ≥ 0, |m|≤ n . Next, from (7), (23) and (40), we have (∂/∂n =n · ∇ ≡ ∂/∂τ on S f , S c ) which yields for every n ≥ 0, |m|≤ n . Similarly, for τ = c and in view of (33) and (21), condition (8) with (24) and (40) leads to c 2n+1 A m n,c + B m n,c − B m n,e = −G m n, > (τ 0 ), for every n ≥ 0, |m|≤ n , while from (9) and (40), we have for every n ≥ 0, |m|≤ n .
In conclusion, the coefficients of the expansions (23), (24), and (40) are given by the solution of the 4 × 4 linear system (42), (44), (45), and (46), i.e., B m n,c = b m n f 2n+1 a n σ n D n − 1 , where D n = (1 + a n )σ n − 1, for every n ≥ 1, |m|≤ n . The special case for n = 0 (m = 0) must be treated separately, since D 0 = 0 (note that a 0 = σ 0 = 0) in the dominators of (47)-(50), providing which concludes the analytical solution of the problem. At this point, we are obliged to make a crucial remark with respect to the undetermined constant coefficient A 0 0, f ∈ R (the rest of the constants are given in terms of A 0 0, f ) that reflects the lack of uniqueness for our problem [15], which is actually due to the fact that Neumann-type conditions are incorporated into the boundary value problem under consideration. As a matter of fact, this is not the case herein, and uniqueness is secured since the interior electric potential (23) is written in terms of the additive constant A 0 0, f , arising from the case n = 0, which vanishes when the Neumann condition applies, while all the other coefficients are evaluated explicitly. However, it is obvious that an additional condition must be imposed in order to calculate this coefficient, which could be the assignment of a specific value for the potential u f at τ = 0 through the general condition (11). For example, the choice u f (0) = 0 fixes the arbitrary additional constant, allowed by the mathematical statement of the problem, i.e., A 0 0, f = 0. Therefore, this peculiarity of the forward EEG problem that appears in the literature is evidently treatable, leading to uniqueness in that sense, and the given solution holds true.

The Single Compartment Limiting Case
In this section, we consider the limiting case where the concentric fluid core vanishes. Then, obviously, the only remaining fields are u c and u e , where we have f → 0 + , hence, for n ≥ 1, {a, a n , D n } → +∞ , while D n = a n , assuming σ f = σ c , therefore σ = 1 and σ n = 1. In this particular case, the coefficient (47) does not correspond to a field; however, it is integrated into the coefficients (48)-(50). Hence, provided that the source functions G m n, > < (τ 0 ) are bounded, relation (52) gives b m n → ±∞ and (47) takes the form for every n ≥ 1, |m|≤ n . As a result, (48) becomes for every n ≥ 1, |m|≤ n . On the other hand, since D n = a n σ n , from (49), it turns out that for every n ≥ 0, |m|≤ n , which is expected as the intermediate cell problem in Ω c has now become an interior one (the corresponding expansion does not contain the term τ −(n+1) for all n ≥ 0). Then, since σ n = 1, from (50), we end up to for every n ≥ 1, |m|≤ n . Finally, for n = 0 and m = 0, due to (54), we obtain which inherits the lack of uniqueness of the Neumann-type boundary value problem and the necessity of additional information, as is readily discussed in the previous section. The result (59) does not cancel either (56) or (58) since they both always hold true, providing zero on each side, due to the fact that G 0 0, > (τ 0 ) = 0 from (39) with (33). To sum up, in this instance, the involved potentials u e and u c are derived from (24) and (40), respectively, i.e., and where the coefficients are given by (56), (58), and (59), while we utilize (39) according to definition (33). The potential fields (60) and (61) coincide with those in the bibliography for the single-layer spherical model, providing an effective analytical validation of our analytical approach.

Numerical Implementation and Results
For the purpose of illustrating the above analytical results and investigating the usefulness of the presented model, we provide a graph (see Figure 2) of the electric potential on the surface of the head, where we consider a dipole Q = (2, 1, 3) × 10 −5 Cm, situated inside the cerebrum Ω c , at a point with coordinates τ 0 = 5 cm, ϑ 0 = π/3 and ϕ 0 = 4π/3. Moreover, we assume that the cerebrum has an outer radius c = 7 cm and electric conductivity σ c = 0.37 S/m [16], while we consider the interior region Ω f with a higher conductivity Indeed, the corresponding Figs. 2(a), 2(b), and 2(c) verify that the presence of an adequately large edema, with different electric conductivity than that of the cerebrum, influences EEG recordings. The differences between cases (a), (b), and (c) are evident in the whole surface of the head, depending on the size of the fluid core, although they seem to be stronger around the peak values of the potential, e.g., along the direction of the current dipole . However, numerical experiments showed that the differences between these three cases become inconsiderable if we assume a value of f σ close to that of c σ , which is expected. Note that the third case (c) almost represents the limiting case analyzed in Section 4, where the fluid core is nonexistent.  As it is demonstrated in Figure 3, the exterior potential e u on a specific point of external surface c S starts to converge after the use of about 15 terms for n in the expansion of the potential (24), i.e., max 15 n = . However, in order to be consistent with high accuracy in practical situations, in Figure 2, we utilized max 25 n = , wherein the use of additional terms does not contribute to the values of the potential.

The Magnetic Field
The magnetic induction field is obtained by employing the Biot-Savart law Due to the spherical symmetry of the present model, in every point of the concentric spherical surfaces S f and S c , the outwardly directed unit normaln coincides with the radial unit vectorτ , i.e.,n =n(τ ) ≡τ , hence, in view of the identity relation (64) yields for every ϑ ∈ [0, π] and ϕ ∈ [0, 2π). We also remind from (35) that where g m for every n ≥ 0 and |m|≤ n , while ν = n and " < ", τ < τ −(n + 1) and " > ", τ > τ , n ≥ 0, depends upon the position τ over the point (τ , ϑ , ϕ ). Moreover, the potentials u f and u c are given by (23) and (40), respectively, which can be rewritten as and Finally, taking the dot product of (67) andτ and provided that the involved surface integrals vanish, therefore we end up tô which is valid for every τ and means that the radial component of B is a function of only the dipole moment Q and its location τ 0 . This is an important result, which has been proved in [9], according to which MEG is independent of any radial variation of the conductivity in a spherical model, which means that it does not depend on the edema region.

The Exterior Magnetic Field
Next, we calculate B in the exterior space Ω e where the total current J t is zero; hence, (18) implies that the magnetic field is both irrotational and solenoidal. Therefore, it accepts the representation for every τ ∈ Ω e . Now, since the magnetic scalar potential U in the exterior space has the asymptotic behavior we obtain hence, from (78) and (79), we arrive at or In this case, where B is also a conservative field (as irrotational), we can choose to integrate along the radial directionτ, from τ to infinity (path independent). Then, the above takes the form which is written as where To conclude, the magnetic induction (79), in view of (85), becomes for every τ ∈ Ω e , hence, for every τ ∈ Ω e . Finally, if we introduce the magnetic potential (85) is written while the magnetic induction (90) takes the form The above results (85) and (90) or (92) and (93) comprise closed-form representations (that do not depend on σ f , σ c ) of the magnetic potential and the magnetic field, respectively, in the exterior of a human head with brain edema, where it is evident that if the point current dipole is radial, the magnetic field vanishes. This is a famous result in MEG research, which can be easily proved with the aid of expression (87). Indeed, a radial dipole assumes the form Q = Qτ, thus considering that the position vector reads τ = ττ, then it is easily shown that Q × τ 0 · τ = Qτ × τ 0 · ττ = Qττ × τ 0 ·τ = 0, since the mixed triple product becomes zero. This outcome is readily substituted to (87), leading to B(τ; τ 0 ) = 0, which verifies the fact that a radial dipole produces no exterior magnetic field.
Obviously, these results for the exterior magnetic field are identical to the corresponding results of the single-compartment spherical model of the cerebrum (without edema), which frequently appear in the ample literature for EEG and MEG.

Conclusions and Discussion
In this work, we provide analytical results for the electric potential and the magnetic field in the interior and the exterior of the head, which is modeled as a non-homogeneous conductor consisting of a spherical fluid core and a concentric spherical shell representing the cerebrum, where the current dipole lies. We observe that the results for the exterior magnetic field are not affected by the existence of the liquid core. However, the interior and exterior electric potential is a function of the fluid core diameter; hence, EEG recordings could be beneficial in confirming the existence of brain swelling. The reduction of the obtained results of the electric potentials for the EEG to those that correspond to the singlelayer problem and the evaluation of the exterior magnetic field that is independent of the interior properties of the cerebrum, provide us with already well-known information, which stands for a valuable criterion of the validity of our method.
The advantage of such analytical solutions and closed-type formulae in comparison with pure numerical methodologies lies in the fact that they readily provide efficient and handy expressions for the evaluation of the implicated fields. Therefore, the validity of numerical solutions can be verified by these analytical or even semi-analytical techniques. On the other hand, fundamental physical laws are derived from analytical methods providing a stable and secure basis for starting a numerical procedure. Consequently, simple analytical methods can co-exist with numerical analysis, cooperating toward the solution of boundary value problems in physical applications of importance.
Our future steps involve further and more sophisticated investigation of the above spherical model in terms of an eccentric spherical edema with respect to the center of the brain tissue and/or a multilayer spherical head model, considering the effect of the skull and the scalp (along with the edema and the brain tissue) only to the EEG problem since the MEG problem does not depend on the skull and scalp conductivities. Finally, future work includes extensive numerical evaluation in order to estimate the effect of the size of the fluid core and the location of the dipole source in EEG recordings, as well as the use of computational experiments with the parameters from real life in order to demonstrate the possibilities of the model to reproduce the situations that are interesting for medicine and biophysics.