Analytic Modeling of Neural Tissue: I. A Spherical Bidomain

Presented here is a model of neural tissue in a conductive medium stimulated by externally injected currents. The tissue is described as a conductively isotropic bidomain, i.e. comprised of intra and extracellular regions that occupy the same space, as well as the membrane that divides them, and the injection currents are described as a pair of source and sink points. The problem is solved in three spatial dimensions and defined in spherical coordinates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(r,\theta,\phi )$\end{document}(r,θ,ϕ). The system of coupled partial differential equations is solved by recasting the problem to be in terms of the membrane and a monodomain, interpreted as a weighted average of the intra and extracellular domains. The membrane and monodomain are defined by the scalar Helmholtz and Laplace equations, respectively, which are both separable in spherical coordinates. Product solutions are thus assumed and given through certain transcendental functions. From these electrical potentials, analytic expressions for current density are derived and from those fields the magnetic flux density is calculated. Numerical examples are considered wherein the interstitial conductivity is varied, as well as the limiting case of the problem simplifying to two dimensions due to azimuthal independence. Finally, future modeling work is discussed.


Introduction
The purpose of this paper is to model the electric potentials in and around a finite volume of excitable tissue that result from externally applied injection current. Our motivation toward quantitative understanding of the distributed electrophysiology of excitable tissue is due to the emergence of magnetic resonance electrical impedance tomography (MREIT) [1]. The contrast in MREIT-as well as in another MR technique, Electrical Properties Tomography (EPT) [2]-depends on the electrical property distribution throughout the region of interest. Briefly, in an MREIT scan, current is injected into an object in concert with the pulse sequence of an MRI scanner. This current will induce a magnetic field [3] whose distribution throughout the entire region can be captured via the phase component of the reconstructed MR images. Electrical conductivity maps may then be constructed from the phase data using the Laplacian of the z component of the induced magnetic field, ∇ 2 B z [4,5]. MREIT has already shown clinical promise, e.g. lesion characterization [6], but it is the possibility of monitoring brain activity with MREIT [7] that especially motivates this study. If MREIT is to be used to detect neural activity it is useful to estimate the influence of MREIT imaging currents on both active and passive tissues. Therefore, we have constructed from first principles an analytic mathematical model of tissue stimulated by injection currents, not unlike that of an MREIT scan.
Excitable tissues are comprised of cells, discrete units through which electric signals may propagate via action potentials [8]. While many have studied and modeled the behavior of individual cells in both sub-and supra-threshold conditions, it is also very important to understand excitability behavior at the tissue level. This approach has been particularly useful in understanding cardiac activity [9]. The bidomain model [10], a generalization of the cable equation [11], has been employed in this area avoiding the discrete constructs of tissue, assuming instead a continuum of two domains, intra-and extracellular, connected by a membrane and that occupy the same volume [12]. Each domain represents an average, then, of all its individual components. MR imaging also necessarily involves averaging over tissues. If we seek to image neural activity using MREIT it is convenient to use a geometrically simple model to predict changes in these images created by neural activity.
Many authors have modeled excitable tissue with the bidomain equations, choosing the coordinate system that most closely resembles the tissue geometry. In circular cylindrical coordinates, Altman and Plonsey modeled a bundle of nerves as an infinite cylinder in an infinite conducting bath, studying first the steady state [13] and transient stimulation [14]. In the former they incrementally increased the realism of their model, going from an isotropic monodomain to an anisotropic bidomain, while in the latter they investigated the effect of fiber diameter on stimulation and impulse propagation. Henriquez et al. [15][16][17] and Trayanova et al. [18] investigated the merits of assuming a single fiber vs. a bundle, i.e. bidomain, of fibers when modeling an infinite cylinder of tissue excited by either a disk or a line source. They showed that the single fiber core conductor model is not an unreasonable approximation of the control region of a large bundle of fibers, but loses its validity toward the periphery of the bundle and is entirely unsatisfactory for small bundles. Plonsey and Barr showed in a two dimensional rectangular framework, except for special cases, the bidomain approach to modeling tissue electrophysiology is not a mere generalization of one dimensional cable theory [19,20]. They found that current flowed very differently in isotropic tissue compared to anisotropic tissue with unequal anisotropy ratios. Roth gave approximate analytic solutions to the problem of bisyncytia with unequal anisotropy ratios [21], using rectangular coordinates. His perturbation method involved expansion in a parameter defined through the anisotropy ratios. He considered two sources: an expanding wave front that was approximated with a step function, and a point source. Trayanova et al. considered the case of bidomain tissue in a uniform electric field, modeling the heart as a sphere of anisotropic tissue with a core of blood [22]. The uniform field meant that they could assume azimuthal independence, leaving only a two dimensional problem in the spherical coordinates r and θ . Heretofore none has studied a three spatial dimension bidomain problem in spherical coordinates.
Our present study is motivated by the need to understand the effect on MREIT images of excitable tissue-specifically, a ganglion excised from the abdomen of a sea slug (Aplysia californica)-affected by injection currents injected through electrodes set into the boundary of its artificial sea water bath [7]. We develop a model that is a dramatic simplification of the actual experiment but which still is novel for its generalization to three spherical dimensions. In and of itself this model will depict basic electrophysiological phenomena and can act as a standard against which numeric simulations such as finite element models (FEM) are held, lending credibility to those in concurrence. Seen in a broader context, this work can serve as the foundation for more and more sophisticated analytic modeling, e.g. nonlinear transmembrane currents and mixed boundary conditions.
In this first study of three dimensional analysis of distributed neural tissue we model the Aplysia abdominal ganglion (AG), known to be electrically coupled by gap junctions [23], as an isotropic bidomain sphere, the artificial sea water bath as an infinite isotropic conducting medium, and the injection currents as source as sink points. We assume isotropic conductivity here for simplicity. However, anisotropy may be the subject of future work, as active tissue is generally anisotropic.

Geometry
Let there be given a sphere of isotropic excitable tissue in a uniform isotropic infinite conducting bath which also contains a point current source and a point current sink. We shall consider this problem in terms of spherical coordinates (r, θ, ϕ) [24]. The sphere of tissue, whose radius is r = a, has its center at the origin. The current source and sink points are distances p + = (r + , θ + , ϕ + ) and p − = (r − , θ − , ϕ − ), respectively, from the origin, as shown in Fig. 1. The current source and sink are in the conducting bath, not in the tissue, i.e. a < r + , r − . The segments R source and R sink are the respective distances from the source and sink to any field point with position vector r = (r, θ, ϕ). The angles γ + and γ − , drawn with a dot-dashed line in Fig. 1 are between p + and r, and p − and r, respectively. Sphere with radius r = a in an infinite conducting medium with current points source and sink at p + = (r + , θ + , ϕ + ) and p − = (r − , θ − , ϕ − ), respectively. The dot-dashed curves labeled as γ are the angles between the points' position vectors and that of a field point r = (r, θ, ϕ)

Bidomain Tissue
The tissue is modeled as a bidomain: two regions-intracellular and extracellularthat occupy the same volume along with the membrane that separates them. Any transmembrane current must be either from the intracellular region to the extracellular region, I m = ∇ · J i , or vice versa, I m = −∇ · J o [10]. From Ohm's law, J = E/ρ, where E is electric field strength and ρ is resistivity [3]. Assuming that E is quasistatic [25] and, further, that there is no tissue capacitance, we may express E in terms of scalar potentials, φ, i.e. E = −∇φ. Thus we arrive at the bidomain equations, where the bidomain potentials, φ i and φ o , are the intra-and extracellular potentials, respectively, and ρ i and ρ o are the corresponding resistivities. Throughout this analysis we assume the membrane to be passive resistor which makes I m depend upon the difference between φ i and φ o , where R m is the membrane resistance times unit area and β is the ratio of the membrane's surface area to volume of the cell. Equations (1a)-(1b) are coupled and must be un-coupled to solve for φ i and φ o by recasting the system in terms of the transmembrane potential, V m , and the monodomain potential, ψ [26], The bidomain potentials are, then, given as To solve for V m , let us subtract Eq. (1b) from Eq. (1a) and then insert Eq. (2), where ρ m = R m /β is the membrane resistance times unit volume. If we define a length constant, λ = √ ρ m /(ρ i + ρ o ), from Eq. (3a) we can immediately see that V m satisfies the scalar Helmholtz equation, To find a relationship that only involves ψ , let us apply the Laplacian operator to Eq. (3b): When we put our expressions from Eqs. (1a)-(1b) for ∇ 2 φ i and ∇ 2 φ o into the right side of Eq. (8), the two terms on that side add to 0, leaving us with the Laplace equation through ψ,

Infinite Medium
External to the tissue the potential, φ e , is given by where φ source and φ sink are the fields due to the current point source and current point sink, respectively, and φ bath is the secondary field [27] which satisfies the Laplace equation,

Transmembrane Potential
The scalar Laplacian of a function f is defined as the divergence of the gradient of f [24] where, in spherical coordinates, u 1 = r, u 2 = θ , u 3 = ϕ, g 1 = 1, g 2 = r 2 , g 3 = r 2 sin 2 (θ ), and √ g = r 2 sin(θ ), whence we arrive at the familiar expression [24] If we apply Eq. (13) to Eq. (7) and assume a product solution of the transmembrane potential V m = R(r)Θ(θ)Φ(ϕ), we can use the method of separation of variables to get three independent, linear, ordinary, second order, differential equations [24], where α = μ(μ + 1). Equation (14a) admits two solutions, i μ and k μ , the modified spherical Bessel functions of the first and second kind of order μ, respectively [28]. We only use i μ , however, because the domain includes the origin, where k μ is singular. The transmembrane potential, then, is where Y ν μ is the tesseral spherical harmonic [29] of degree μ and order ν, and a μν is the coefficient determined from the boundary conditions.

External Potential
The potential from the external bath also satisfies the Laplace equation, but its domain does not include the origin, so we may immediately write its solution as where c μν is determined by the boundary conditions. The potentials due to the current source and sink points with magnitude I o are given as [10] where ρ e is the resistivity of the conducting bath, and R source and R sink are the distances from their respective points to any point (r, θ, ϕ) in the problem domain, shown in Fig. 1. To satisfy the boundary conditions, φ source and φ sink must be written in terms of r, θ , and ϕ, the derivation of which can be found in the Appendix.

Boundary Conditions
We have three boundary conditions [30] at the tissue-bath interface, where r = a, through which we will determine the unknown coefficients, a μν , b μν , and c μν . They are continuity of external and extracellular potentials, continuity of normal current between bath and interstitium, and no intracellular normal current whose solutions yield where and completely determining all potential fields.

Current Densities
Current density is proportional to the negative gradient of the scalar electric potential, J = −ρ −1 ∇φ. The gradient of a function f is given as [31] In spherical coordinates u i and g i are the same as in Eq. (12) and the unit vectors are a 1 = r, a 2 = θ , and a 3 = ϕ which gives us Applying Eq. (27) to Eqs. (4a)-(4b) and (10) and dividing by the resistivities of their respective domains gives us expressions for the current densities throughout our problem. They are J e (r, θ, ϕ) = c μν ρ e r 2+μ + where and Here g < +,− = min(r, r +,− ), g > +,− = max(r, r +,− ), and Γ is the gamma function [32].

Magnetic Flux Density
Here we summarize the theory behind our numeric calculation of the magnetic flux density, B, due to J. Within a current-carrying volume, B may be calculated from the Biot-Savart law [33], thus every point in B requires integration over the entire current-carrying volume.
Since the numerator of the integrand of Eq. (34) is a cross product the individual components of B in cartesian coordinates are given as In an MREIT scan we are only concerned with the component of the magnetic field that is along the axis of the bore of the MRI scanner, i.e. B z ; so, we will focus on that now but the following would apply to B x and B y as well. We can see that B z is the difference of two integrals, If we define the x and y components of a function G as then we can define B z as a difference of two convolutions [34], Convolution becomes multiplication in the Fourier domain [34] thus B z finally is calculated as where F is the Fourier transform [34].

Numeric Examples and Discussion
We now demonstrate our modeling with some simple examples in which we show the effect of varying ρ o . Let us assume the tissue radius is a = 2 mm, the point current source is at p + = (5, π 2 , 0), and the point current sink is at p − = (5, π, 0). Plainly said, there is a cathode 2 mm to the right of, and an anode 2 mm directly below, a ball of tissue in an otherwise uniform ocean of artificial sea water. We choose resistivities in a range known to be biologically realistic [10], specifically ρ e = 0.29 Ωm, ρ i = 0.19 Ωm, and ρ m = 0.15 Ωm. The source and sink are positive and negative, respectively, with equal magnitude I 0 = 1 mA. The upper bound was chosen that the solution stabilized to 5 decimal places, i.e. μ = 10. All of these inputs are summarized in Table 1     In the bottom graphs we show the magnetic flux density to which J e , J i , and J o give rise. We composed a program in MATLAB (The MathWorks, Inc., Natick, MA) that employs the fast Fourier and inverse Fourier transforms to compute B z from the simulated J data throughout a 3 mm × 3 mm × 3 mm cube whose origin is the same as the sphere's. The J sampling was 13 slices, evenly spaced along the z axis, each containing 128 × 128 data points. All of these results are as expected. But for the dashed white line, the sphere of tissue is indistinguishable from the surrounding bath in Fig. 2. When the interstitium has one tenth the bath's resistivity (Fig. 3) the current can be seen to go toward the sphere resulting in a large B z at the sphere's edge in the fourth quadrant of the field of view. Accordingly, when the interstitium is ten times as resistive as the bath (Fig. 4), the current flows mostly around it. In terms of the B z field, there appears a faint glow around the sphere's edge, surrounding a region of relative darkness, due to this flow pattern. In both cases of ρ o = ρ e the tissue is clearly visible in all three field types, φ, J, and B z . As a limiting example, let us move the source point to be directly above the tissue such that p + = (5, 0, 0). This is an axially symmetric problem, identical to our earlier work [35], and we get the same results, shown in Figs In this article we set out to model the electromagnetic fields in and around a volume of neural tissue stimulated by current that is injected in close proximity to it. The geometry selected is expected in an in vitro MREIT scan, where an AG may be submerged in a bath of artificial seawater contained in a cylindrical sample chamber that has injection current ports on opposing sides [7]. We note that the effect of the applied external field is to simultaneously depolarize and hyperpolarize portions of the simulated tissue nearby the current sources. If a portion of tissue is sufficiently depolarized to form an action potential it may propagate throughout the tissue from these regions. It has been suggested [36] that modest depolarizations or hyperpolarizations caused by weak external currents applied to the skull are sufficient to excite or inhibit neural excitability in brain structures. More complex models of the tissue and field geometry used here may prove useful methods of exploring the mechanisms of such neuromodulation techniques. In these numeric examples we have held the source and sink points to be equidistant from the sphere of tissue, r + = r − . This gives the problem a symmetry about the lines y = −x in Figs. 2, 3, and 4 and y = 0 in Fig. 5. However, the AG is smaller than the diameter of the sample chamber; so, it will not necessarily be directly between the ports, spoiling this axial symmetry. Thus a complete three dimensional treatment of this type of problem is finally required.
From Eqs. (4a)-(4b) we can see that each region, intracellular and extracellular, has a monodomain component ψ that obeys the Laplace equation. Plots of this potential produce results similar to those of Rush and Driscoll [37,38]. They solved for the electric potential in a brain from electrodes placed directly on a scalp, modeling the brain, skull, and scalp as different layers of monodomain tissues, i.e. a sphere encased in a thin shell of bone which was itself encased in a thin shell of skin. We could amend our model to include similar surrounding layers, each with its own expressions for φ and J and coefficients determined from the boundary conditions. The boundary conditions themselves would change, e.g. the interstitium would have continuity of potential and normal current with the skull rather than with artificial sea water. Such changes would be appropriate for a model on the scale of e.g. a dog's head [39].
We have modeled both domains as being ohmic, i.e. their impedivities z = ρ are only real valued, but z can be made complex by introducing a frequency dependence [34]. In their extensive literature review [40] and experimental measurements [41], Gabriel et al. show that most tissues have frequency dependent electrical properties. More recently, Bédard et al. [42] and Bazhenov et al. [43] explored frequency dependence in local field potentials. In a series of theoretical papers Bédard and Destexhe provide a general framework for modeling electromagnetic fields in brain tissue without assuming the interstitium to be purely resistive. Absent those assumptions, they developed a generalized formalism of current source density analysis with the goal of relating the extracellular potential to current sources in the tissue [44], considering monopolar sources, dipolar sources, and combinations thereof. Next they incorporated frequency dependent extracellular and intracellular impedivities, z o and z i , to generalize the cable theory [11] for neurons embedded in a complex interstitium [45]. They showed that z o and z i have a non-trivial impact on the properties of neurons, e.g. voltage attenuation with distance and the spectral profile of V m . Finally, they calculated the magnetic fields generated from a current-carrying neuron and, using superposition, a population of neurons [46]. They showed that since the electrical properties of neural tissue impact the transmembrane and axial currents of a neuron, they will necessarily also impact the magnetic fields these currents create. By contrast, in our study we have concerned ourselves with the interaction of neural tissue and an aphysiologic stimulus. The effects of this stimulus will naturally also depend upon complex tissue properties, but over a larger scale determined by the stimulus geometry. Future work should explore the impact non-ohmic impedivities have on a tissue interactions with applied external fields.
Let us now consider possible next steps to build on this first study. As well as modeling tissue properties as complex, it should also be possible to examine the transient behavior of excitable tissue using a Hodgkin-Huxley-like model [47]. These analytical models of spheres could then be validated and used to estimate the scale of changes expected in MREIT images due to different neural activity patterns. Fig. 7 Geometry of the point current source at p + = (r + , θ + , ϕ + ) in relation to the origin (0, 0, 0) and a field point r = (r, θ, ϕ), with R source being the distance between r and p + and whose inverse has the form of the generating function of the Legendre polynomial [32], allowing us to write 1 r 2 + r 2 + − 2rr + cos(γ + ) where P μ is the Legendre polynomial of order μ. The spherical law of cosines [29] lets us write γ in terms of θ and ϕ, cos(γ + ) = cos(θ ) cos(θ + ) + sin(θ ) sin(θ + ) cos(ϕ − ϕ + ), which, from spherical harmonic addition theorem [48], allows us to expand P μ in terms of Y ν μ , The analysis for a current point sink proceeds analogously. From Eqs. (19a)-(19b), the potentials for current points source and sink finally are given in spherical coordinates as φ source (r, θ, ϕ) = I 0 ρ e ∞ μ=0 μ ν=−μ (g < + ) μ (g > + ) μ+1 φ sink (r, θ, ϕ) = −I 0 ρ e ∞ μ=0 μ ν=−μ (g < − ) μ (g > − ) μ+1