Finite element analysis of neuronal electric fields: the effect of heterogeneous resistivity
Simulation of extracellular fields is one of the substantial methods used in the area of computational neuroscience. Its most common usage is validation of experimental methods as EEG and extracellular spike recordings or modeling of physiological phenomena which can not be easily determined empirically. Continuous experimental work has been re-raising the importance of polarization effects between neuronal structures to neuronal communication. As this effects rely on very small potential changes, better modeling methods are necessary to quantify the weak electrical fields in the microscopic scale in a more realistic way.
An important factor of influence on local field effects in the hippocampal formation is the heterogeneous resistivity of extracellular tissue. The vast majority of modeling studies consider the extracellular space to be homogeneous while experimentally, it has been shown that the stratum pyramidale has two times higher resistivity than other hippocampal layers. Common simulation methods for extracellular electrical fields based on the point source approximation are bound to describe the resistance of the space with a single scalar. We propose that models should be based on the space- and time-dependent Maxwell equations (Partial Differential Equations, PDEs) in order to account for heterogeneous properties of the extracellular space and specific arrangements of neurons in dense hippocampal layers.
To demonstrate the influence of heterogeneous extracellular resistivity and neuronal spatial orientation on modeling results, we combine solutions of classical compartment models with spatiotemporal PDEs solved by the Finite Element Method (FEM). With the help of these methods, we show that the inclusion of heterogeneous resistivity has a substantial impact on voltages in close proximity to emitting hippocampal neurons, substantially increasing the change in extracellular potentials compared to the homogeneous variant.
Numerous computational studies have investigated time-varying currents in homogeneous extracellular space [1, 2, 3] as well as the role of neuronal morphology in uniform electric field stimulation [4, 5, 6]. However, most studies use models of the extracellular milieu that may not be accurately applied when modeling brain structures with highly diverse extracellular resistivity and neuronal arrangement.
As emphasized by Lopez-Aguado and Bokil [7, 8] it has been traditionally neglected that currents propagate in all directions in an extracellular medium and that inward and outward currents originate from tissue regions having large resistivity differences. A usual argument for this approach is that resistivity influences extracellular electrical fields minimally and that extracellular space can be assumed as homogeneous in current-source density (CSD). However, extracellular non-homogeneous resistivity has been shown experimentally in several regions of the brain, for example, in the hippocampus and the cerebellum . Additionally, tissue swelling has been observed after intense neural activity that in turn could lead to an increase in extracellular resistivity . Previous studies [1, 11, 12, 13] have attempted to examine how extracellular electrical fields affects neuronal activity although with the help of quasi-static approximation and an assumed homogeneous tissue resistivity.
To quantify the effect of extracellular heterogeneous resistivity and neuronal spatial distribution on strength of neural fields, we are proposing a simple modeling pathway to couple compartment-based neural models with the COMSOL Multiphysics simulation environment. Here, we solve time-dependent Maxwell’s equations using the FEM to analyze the change of electrical fields as it occurs in the extracellular space surrounding neurons.
Our results indicate that inhomogeneous resistivity of the extracellular milieu significantly influences the change of extracellular potentials (EPs) in the hippocampus. By computing the resulting voltage change due to outgoing transmembrane currents of an exemplar CA1 pyramidal cell model in homogeneous and inhomogeneous extracellular resistivity, we observed a maximal difference in EP change of 60% in the hippocampal pyramidal layer. Furthermore, the here proposed method offers the possibility to efficiently simulate the effects of superimposed extracellular potentials created by neurons in divergent positions relative to each other. We will also discuss advantages and drawbacks of the used method and propose alternatives and possible improvements towards more realistic modeling of electrical fields of the brain.
Ii Materials and Methods
Ii-a Pyramidal neuron model
To demonstrate the effect of trans-membrane currents on the effect of extracellular fields, a Hodgkin-Huxley (HH)-like hippocampal CA1 pyramidal cells was adapted from  with the addition of a current from  and a current from .
Our model contained one somatic (r = 10 ), 20 dendritic (r = 1–3 , l = 5 ), two axon initial segments, (r = 0.5 , l = 60 ), (r = 0.5 , l = 60 ), and finally 32 axonal compartments (r = 0.5 , l = 5 ), see Figure 1. Note that here refers to the part of AIS from 0 to 60 , and from 60 to 120 from the cell body. Additionally, compartments contain high channel density . 20 dendritic compartments are used for the dendritic tree, while additional compartments are included in the branching analysis. The channel density was varied between 986 and 2943 (corresponding to ). The following conductance values were used: = 8, = 5, = 1 for the soma; = 10, = 15, = 0.8, = 0.625, = 0.07, = 0.2 for the dendrite; = 50, = 10, = 10 for the and = 9, = 10, = 10 for the rest of the axon compartments. The leakage conductance was set to 0.1 for all of the compartments. Finally, the equilibrium potentials were set to = 60 mV, = -85 mV, = -43 mV and = -65 mV.
Ii-B Creation of morphology and model coupling
The three-dimensional neuronal geometry was constructed in COMSOL Multiphysics 4.3 with the help of the interface to MATLAB (“LiveLink”) by morphological additions and boolean unification of simple geometric volumes. As we aim to represent the 3D-morphology as an exact counterpart of the compartmental model, each section is recreated as a cylinder with the same length and diameter as in the compartmental definition. If the cylinders were added on top of each other with no change of the rotation vector, as shown in Figure 1, no joining geometrical primitives were added in between them. If the rotation differs, as for example in Figure 4, a sphere is added in between the cylinders, followed by a removal of the interior boundaries of both the cylinders and the sphere. The reason for constructing the geometry in this way is that the mesh engine otherwise respects the internal boundaries such that the resulting mesh becomes unnecessary complex. In some cases the occurrence of multiple internal boundaries can even make the meshing procedure abort.
Ionic currents of single neuronal compartments model were determined by solving Hodgkin-Huxley Ordinary Differential Equations (HH ODEs) specified in equations 1 to 6 using a Runge-Kutta algorithm from the MATLAB ODE-suite (Mathworks) with a constant step-size . The sum of all currents per compartment for each time step of the simulation is afterwards normalized to an absolute value of units and stored in a matrix of size .
Next, by using the interface to COMSOL, each row of the matrix was mapped adequately to the corresponding cylindrical domain as a boundary current source . Note that we hereby assume that transmembrane currents are the only cause of change of extracellular potential, which is surely not the case in a real neuron, as for example synaptic calcium-mediated currents are suspected to contribute to a large fraction of the extracellular signature .
The extracellular volume was modeled by cylindrical objects covering the neuronal morphology with either homogeneous (0.3 S/m) or heterogeneous resistivity as shown in Figure 2A. In this case we constructed the extracellular volume by taking the union of cylindrical objects with increasing resistivity in the -axis, according to . It is assumed that the conductivity of extracellular tissue is frequency-independent in the used range of neural activity (10–100 Hz) .
Ii-C Electrostatic formulation and the Finite Element Method
To simulate non-homogeneous distributions of electrical fields produced by single neurons and neuronal networks, we used an electrostatic formulation of Maxwell’s equations discretized with finite elements. Chiefly, in the finite element approach, Maxwell’s equations are solved by discretizing the incorporated volumes (in this case the neuronal compartments) into finite tetrahedral volume elements . We seek the electric field intensity in terms of the electric scalar potential ,
The relevant dynamic form of the continuity equation with current sources is given by
with and the current density and electric charge density, respectively. Further constitutive relations include
and Ohm’s law
in which denotes the electric flux density. Finally, Gauss’ law states that
Rewriting the electric charge density using Gauss’ law together with the constitutive relation (9) and finally applying the gauge condition (7) twice we arrive at the time-dependent potential formulation
This is the formulation used in COMSOL Multiphysics . The values for the electric conductance and the relative permittivity were obtained from . The source currents were computed from the compartmental model as described above.
The formulation (13) is efficiently solved by COMSOL’s “Time discrete solver”, which is based on the observation that the variable satisfies a simple ODE. Solving for in an independent manner up to time , it is then straightforward to solve a single static PDE to arrive at the potential itself.
As for boundary conditions we took homogeneous Neumann conditions (electric isolation) everywhere except for in a single point which we choosed to be ground (). In all our simulations this point was placed at the axis of rotation of the enclosing cylindrical extracellular space, and underneath the neuronal geometry. This procedure ensures that the formulation has a unique solution (it is otherwise only specified up to a constant).
A tetrahedral mesh was applied to discretize space (using the “finest” mesh setting; resolution of curvature 0.2, resolution of narrow regions 0.8). The simulations were verified against coarser mesh settings in order to ensure a practically converged solution. As a final note, in the Time discrete solver the time step was set to the same step size as used in the ODE-based solution of the Hodgkin-Huxley equations, thus ensuring a correct transition between both simulation environments.
Iii-a The effect of heterogeneous extracellular fields
In order to examine the influence of heterogeneous extracellular space, we first constructed a three-dimensional active neuron model in an homogeneous and heterogeneous extracellular milieu (Figure 2). Neurons in the hippocampus have an intricate spatial orientation that propitiates strong field potentials: high density of pyramidal cell dendrites running in parallel in the stratum radiatum (SR), densely packed pyramidal cell somas in the stratum pyramidale (SP) while pyramidal cell axons run almost in parallel or crossing each other in SP or stratum oriens (SO).
We measured the voltage on four defined point probes placed parallel to the dendrites, soma, and axon terminal compartments during the peak of an action potential (AP) of 80mV, by varying the distance between active neuron and point probes from 1 to 80 (Figure 2A).
In the first set of simulations, we analyzed the effect of the aforementioned four neuronal regions on the defined point probes assuming a widely accepted homogeneous resistivity of 350 cm [2, 3]. Note that here applied boundary currents of three-dimensional neuronal compartments correspond to the peak transmembrane current during an AP. By doing so, the peak voltage of 0.25 mV was obtained in the point probe parallel to the compartment, followed by soma, dendrite and axon terminal. Consequently, the heterogeneous case was examined by placing the neuron and point probes into a heterogeneous extracellular space representing hippocampal spatial order, where resistivity values for different strata were obtained from  and shown in Figure 2A. The active neuron was positioned in the center of SP and point probes were moved along the x axis of the extracellular space. In the case of non-homogeneous extracellular resistivity (Figure 2C) the largest voltage change was measured in point probes parallel to , although the values in point probes placed in SP were 60% larger in close distance than in the homogeneous extracellular space scenario and 28% higher considering the spatial mean over the total distance. The point probes placed in SO and SR were as well affected by the higher restivity of the pyramidal layer, showing an average increase of 4% in parallel to the axon and 7% in parallel to the dendrite. The appreciable difference between the voltage changes suggests that non-homogeneous distribution of resistivity is an important aspect of extracellular field effects.
Iii-B Analysis in respect to spatial orientation and neural morphology
According to the superposition principle electrical fields (EFs) are produced by summation of single neuronal activity. Thus, the mutual interaction between neuron and field is strongly modulated by the spatial orientation of neuronal assemblies . To analyze EFs considering a realistic spatial order, we began simulating two neurons with parallel axons (Figure 3A) where represents an active structure with boundary current source mapped to its surface while is a passive measurement structure. We measured at four points along the axon of and detected the maximal voltage amplitude of 0.28 mV. Note that here corresponds to extracellular voltage measured on the cell membrane. Next we tested the influence of four neighboring cells with parallel axons on the (Figure 3B). In this case, having four active neuronal neighbors (distance between the axons =10 ) firing non-synchronously, a maximal voltage of 0.49 mV in homogeneous and 0.68 mV in heterogeneous extracellular medium of was computed. Synchronous firing of adjacent neurons produced a peak voltage of 0.62 mV in homogeneous and 0.84 mV in heterogeneous medium of . When axons crossed each other at the region (Figure 3C), we observed a maximal voltage amplitude of 0.89 mV in the of during asynchronous activity of the neighboring cells, whereas synchronous activity produced 1.41 mV in heterogeneous medium. These results suggest that for analyzed neuronal arrangement soma and have considerable influence on the strength of neural fields for small distances ( 10 ) and that values in heterogeneous extracellular space are in mean 25% larger than in homogeneous extracellular medium when averaged over time.
Furthermore, we simulated the influence of pyramidal neuron position within different hippocampal layers on the strength of neural fields (Figure 4). We positioned somas in SP with the smallest distance between the two axon initial segments confined to SO. This geometrical arrangement was motivated by confocal images of two proximal pyramidal neurons located parallel to each other with an axonal bend starting about 50 behind the soma (Figure 4A). As known from previous analytical studies  the bend of axonal structures generally amplifies endogenous fields around neurons. We were able to confirm this effect in our simulation as we simulated the field effect resulting from trans-membrane currents of an active neuron measured on the parallel neuron. In this scenario, peak voltage of 0.35 mV was calculated at the position of while the smallest voltages were registered at the passive neurons between the dendrites (0.12 mV) in SR and axon terminals (0.16 mV). If were confined to SP, the potential computed at the passive neuron was 25% larger than in the case in which were positioned in the SO (Figure 4B).
In this work, we use the Finite Element Method coupled to the HH-equations to simulate how neuronal geometry, arrangement and heterogeneous extracellular properties affect the strengths of neural fields. We first show that there is a difference in voltage transients produced by firing of neighboring cells if the extracellular space is considered to be homogeneous or heterogeneous. Additionally, we demonstrate that the spatial orientation of specific cellular compartments is an important determinant of the strength of neural fields.
In our computations, the highest change of extracellular potential arises in the pyramidal layer, in proximity to the compartment. Action potentials are generally initiated in the AIS due to a higher density of voltage-gated channels [17, 24, 25], which is reflected by the parameter settings for this compartment in the occupied neuronal model. Additionally, the spatial arrangement of hippocampal pyramidal neurons also propitiates the proximity of in both SP and SO (see Figure 4, ).
Holt and Koch  showed that interactions near cell bodies are more important than interactions between axons by using standard one-dimensional cable theory and volume conductor theory. Another study  reported a 4.5 mV change (in the AIS) caused by extracellular interactions, a change more than 4 larger than what we have found in our simulations of heterogeneous extracellular media. Additionally, by applying analytical methods, Bokil et al  have shown that, in the olfactory system, an action potential of 100 mV amplitude in one axon could produce depolarization in other axons in the bundle sufficient to initiate an action potential. However, using FEM simulations it was not possible to trigger spikes in neurons solely by electrical fields mediated by in- and outflow of trans-membrane currents, in agreement with the work by Traub and colleagues (in which maximal voltages in a sink axon during synchronized activity of four neighboring neurons is ).
Hence, simulations relying on the point source approximation to describe hippocampal neural fields may be distorted as the potential change is more than 28% greater in stratum pyramidale and in average 6% greater in other hippocampal layers if a heterogeneous tissue is used instead of a homogeneous one [7, 1].
For varying extracellular resistivity a numerical procedure referred to as ‘the method of images’  has been proposed as an extension to the point source equation. Although this method has been used in computations of extracellular action potentials in a previous study , its practical use is limited to a rather low number of resistivity layers and non-complex geometry.
The requirement of FEM to model passive current flow between neurons was suggested elsewhere [2, 3], but implementation issues and the lack of adequate software tools may have precluded its usage in the past.
Extracellular resistivity could also contribute to the amplitude of local field potentials (LFP), and in fact, in vitro LFP registered in hippocampal slices are greater in the SP than in other layers . Interestingly, we observed that hippocampal slices in interface-type recording chambers (where slices are not completely submerged in the buffer solution [31, 32]) are more than 60% less conductive than in the chambers where slices are submerged (unpublished results). This could help explaining why LFP recorded in interface-type cambers are far greater than in submerged-type counter parts [16, 32].
Here we show that the use of FEM software (COMSOL Multyphysics), with an interface to Hodgkin and Huxley-type ODE model in Matlab, is a powerful tool to verify macroscopic effects of neural fields. However, this approach may fail to simulate small nuances of extracellular interactions (e.g. ion channels apposing active compartments may suffer more from the effect of passive current flow than channels in the opposite side). Nonetheless, simulations using PDEs could increase the level of realistic models of membrane dynamics. The idea of translating membrane dynamics to PDE was proposed by Hodgkin and Huxley themselves , and later pursued by Kashef and Bellman . However, implementing membrane dynamics in combination with the Maxwell equation interface of the PDE simulator has been, to our experience, quite cumbersome. For example, due to the relatively high computational cost of the solution phase of the finite element method, a representation of fully reconstructed morphologies with volumetric cylindrical elements was not possible. The high amount of (small) cylindrical elements required to accurately model complex dendritic branches typically caused the meshing algorithm to break.
The software used in our simulations (COMSOL Multiphysics) has helped to popularize the use of FEM in neuroscientific problems [35, 36]. However, the software is geared towards industrial applications and the steep learning curve associated to the adaptation of COMSOL to neuroscience related problems may preclude its widespread usage by the neuroscience community. Possible combination of Neuron Simulation Environment (NSE) with a FEM simulator would offer great possibilities by allowing both realistic neural morphologies and extracellular space modeling. However, currently the existence of a direct interface between NSE and a FEM simulator has not been reported; future studies should focus on this task. Furthermore, hippocampal extracellular stimulation has been proposed as therapeutic strategy to control several disorders, including temporal lobe epilepsy [37, 38]. Developing models that consider as many as possible physical tissue properties would help to improve extracellular stimulation in epilepsy patients. In summary, our work adds to the recently published studies attempting to reveal important parameters determining the strength of extracellular electric fields. Models that do not use space and time-dependent differential equations when modeling neuronal interactions may have failed to replicate the changes in measured voltages caused by passive current flow in heterogeneous extracellular tissue, especially when more than two neurons are modeled simultaneously.
This work was supported by Kjell and Märta Beijers Foundation and the BMWF Marietta Blau stipend. PB and SE were supported by the Swedish Research Council within the UPMARC Linnaeus center of Excellence.
-  C. C. McIntyre and W. M. Grill, “Excitation of central nervous system neurons by nonuniform electric fields,” Biophysical Journal, vol. 76, no. 2, pp. 878–888, Feb. 1999, PMID: 9929489. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/9929489
-  C. Bédard, H. Kröger, and A. Destexhe, “Modeling extracellular field potentials and the frequency-filtering properties of extracellular space,” Biophysical Journal, vol. 86, no. 3, pp. 1829–1842, Mar. 2004, PMID: 14990509. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/14990509
-  ——, “Model of low-pass filtering of local field potentials in brain tissue,” Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics, vol. 73, no. 5 Pt 1, p. 051911, May 2006, PMID: 16802971. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/16802971
-  C. Y. Chan, J. Hounsgaard, and C. Nicholson, “Effects of electric fields on transmembrane potential and excitability of turtle cerebellar purkinje cells in vitro,” The Journal of Physiology, vol. 402, pp. 751–771, Aug. 1988, PMID: 3236254. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/3236254
-  T. Radman, Y. Su, J. H. An, L. C. Parra, and M. Bikson, “Spike timing amplifies the effect of electric fields on neurons: implications for endogenous field effects,” The Journal of Neuroscience: The Official Journal of the Society for Neuroscience, vol. 27, no. 11, pp. 3030–3036, Mar. 2007, PMID: 17360926. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/17360926
-  T. Radman, R. L. Ramos, J. C. Brumberg, and M. Bikson, “Role of cortical cell type and morphology in subthreshold and suprathreshold uniform electric field stimulation in vitro,” Brain Stimulation, vol. 2, no. 4, pp. 215–228, 228.e1–3, Oct. 2009, PMID: 20161507. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/20161507
-  L. López-Aguado, J. M. Ibarz, and O. Herreras, “Activity-dependent changes of tissue resistivity in the CA1 region in vivo are layer-specific: modulation of evoked potentials,” Neuroscience, vol. 108, no. 2, pp. 249–262, 2001, PMID: 11734358. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/11734358
-  H. Bokil, N. Laaris, K. Blinder, M. Ennis, and A. Keller, “Ephaptic interactions in the mammalian olfactory system,” The Journal of Neuroscience: The Official Journal of the Society for Neuroscience, vol. 21, no. 20, p. RC173, Oct. 2001, PMID: 11588203. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/11588203
-  Y. C. Okada, J. C. Huang, M. E. Rice, D. Tranchina, and C. Nicholson, “Origin of the apparent tissue conductivity in the molecular and granular layers of the in vitro turtle cerebellum and the interpretation of current source-density analysis,” Journal of Neurophysiology, vol. 72, no. 2, pp. 742–753, Aug. 1994, PMID: 7983532. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/7983532
-  A. M. Autere, K. Lamsa, K. Kaila, and T. Taira, “Synaptic activation of GABAA receptors induces neuronal uptake of ca2+ in adult rat hippocampal slices,” Journal of Neurophysiology, vol. 81, no. 2, pp. 811–816, Feb. 1999, PMID: 10036281. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/10036281
-  C. Gold, D. A. Henze, C. Koch, and G. Buzsáki, “On the origin of the extracellular action potential waveform: A modeling study,” Journal of Neurophysiology, vol. 95, no. 5, pp. 3113–3128, May 2006, PMID: 16467426. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/16467426
-  C. Gold, D. A. Henze, and C. Koch, “Using extracellular action potential recordings to constrain compartmental models,” Journal of Computational Neuroscience, vol. 23, no. 1, pp. 39–58, Aug. 2007, PMID: 17273940. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/17273940
-  C. A. Anastassiou, S. M. Montgomery, M. Barahona, G. Buzsáki, and C. Koch, “The effect of spatially inhomogeneous extracellular electric fields on neurons,” The Journal of Neuroscience: The Official Journal of the Society for Neuroscience, vol. 30, no. 5, pp. 1925–1936, Feb. 2010, PMID: 20130201. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/20130201
-  P. F. Pinsky and J. Rinzel, “Intrinsic and network rhythmogenesis in a reduced traub model for CA3 neurons,” Journal of Computational Neuroscience, vol. 1, no. 1-2, pp. 39–60, Jun. 1994, PMID: 8792224. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/8792224
-  R. N. Leao, K. Svahn, A. Berntson, and B. Walmsley, “Hyperpolarization-activated (I) currents in auditory brainstem neurons of normal and congenitally deaf mice,” The European Journal of Neuroscience, vol. 22, no. 1, pp. 147–157, Jul. 2005, PMID: 16029204. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/16029204
-  R. N. Leão, H. M. Tan, and A. Fisahn, “Kv7/KCNQ channels control action potential phasing of pyramidal neurons during hippocampal gamma oscillations in vitro,” The Journal of Neuroscience: The Official Journal of the Society for Neuroscience, vol. 29, no. 42, pp. 13 353–13 364, Oct. 2009, PMID: 19846723. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/19846723
-  L. M. Palmer and G. J. Stuart, “Site of action potential initiation in layer 5 pyramidal neurons,” The Journal of Neuroscience: The Official Journal of the Society for Neuroscience, vol. 26, no. 6, pp. 1854–1863, Feb. 2006, PMID: 16467534. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/16467534
-  K. C. Buszaki G, Anastassiou CA, “Field effects in the CNS play functional roles,” Nature Neurosciecne, vol. 6, 2012. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/22595786
-  H. Lindén, K. H. Pettersen, and G. T. Einevoll, “Intrinsic dendritic filtering gives low-pass power spectra of local field potentials,” Journal of Computational Neuroscience, vol. 29, no. 3, pp. 423–444, Dec. 2010, PMID: 20502952. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/20502952
-  J. Jin, The Finite Element Method in Electromagnetics, 2nd ed. New York: John Wiley & Sons, 2002.
-  AC/DC Module User’s Guide, Comsol, 2012, version 4.3.
-  F. Fröhlich and D. A. McCormick, “Endogenous electric fields may guide neocortical network activity,” Neuron, vol. 67, no. 1, pp. 129–143, Jul. 2010, PMID: 20624597. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/20624597
-  D. Tranchina and C. Nicholson, “A model for the polarization of neurons by extrinsically applied electric fields,” Biophysical Journal, vol. 50, no. 6, pp. 1139–1156, Dec. 1986, PMID: 3801574. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/3801574
-  B. P. Bean, “The action potential in mammalian central neurons,” Nature Reviews. Neuroscience, vol. 8, no. 6, pp. 451–465, Jun. 2007, PMID: 17514198. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/17514198
-  M. H. P. Kole and G. J. Stuart, “Is action potential threshold lowest in the axon?” Nature Neuroscience, vol. 11, no. 11, pp. 1253–1255, Nov. 2008, PMID: 18836442. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/18836442
-  S. Ramón y Cajal, Histologie du Systeme Nerveux de l’Homme et des Vertebretes, 1911.
-  G. R. Holt and C. Koch, “Electrical interactions via the extracellular potential near cell bodies,” Journal of Computational Neuroscience, vol. 6, no. 2, pp. 169–184, Apr. 1999, PMID: 10333161. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/10333161
-  R. D. Traub, F. E. Dudek, C. P. Taylor, and W. D. Knowles, “Simulation of hippocampal afterdischarges synchronized by electrical interactions,” Neuroscience, vol. 14, no. 4, pp. 1033–1038, Apr. 1985, PMID: 2987752. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/2987752
-  E. Weber, Electro Magnetic Fields: Theory, Applications. New York: John Wiley & Sons, 1950.
-  A. Fisahn, F. G. Pike, E. H. Buhl, and O. Paulsen, “Cholinergic induction of network oscillations at 40 hz in the hippocampus in vitro,” Nature, vol. 394, no. 6689, pp. 186–189, Jul. 1998, PMID: 9671302. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/9671302
-  Y. Inaba and M. Avoli, “Volume-conducted epileptiform events between adjacent necortical slices in an interface tissue chamber,” Journal of Neuroscience Methods, vol. 151, no. 2, pp. 287–290, Mar. 2006, PMID: 16143402. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/16143402
-  N. Hájos, T. J. Ellender, R. Zemankovics, E. O. Mann, R. Exley, S. J. Cragg, T. F. Freund, and O. Paulsen, “Maintaining network activity in submerged hippocampal slices: importance of oxygen supply,” The European Journal of Neuroscience, vol. 29, no. 2, pp. 319–327, Jan. 2009, PMID: 19200237. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/19200237
-  A. L. HODGKIN and A. F. HUXLEY, “A quantitative description of membrane current and its application to conduction and excitation in nerve,” The Journal of Physiology, vol. 117, no. 4, pp. 500–544, Aug. 1952, PMID: 12991237. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/12991237
-  B. Kashef and R. Bellman, “Solution of the partial differential equation of the Hodgkin-Huxley model using differential quadrature,” Mathematical Biosciences, vol. 19, no. 1-2, pp. 1–8, 1974.
-  T. C. Zhang and W. M. Grill, “Modeling deep brain stimulation: point source approximation versus realistic representation of the electrode,” Journal of Neural Engineering, vol. 7, no. 6, p. 066009, Dec. 2010, PMID: 21084730. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/21084730
-  N. Yousif, N. Purswani, R. Bayford, D. Nandi, P. Bain, and X. Liu, “Evaluating the impact of the deep brain stimulation induced electric field on subthalamic neurons: a computational modelling study,” Journal of Neuroscience Methods, vol. 188, no. 1, pp. 105–112, Apr. 2010, PMID: 20116398. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/20116398
-  P. H. Stypulkowski, J. E. Giftakis, and T. M. Billstrom, “Development of a large animal model for investigation of deep brain stimulation for epilepsy,” Stereotactic and functional neurosurgery, vol. 89, no. 2, pp. 111–122, Apr. 2011, PMID: 21336007.
-  H. Luna-Munguía, A. Meneses, F. Peña-Ortega, A. Gaona, and L. Rocha, “Effects of hippocampal high-frequency electrical stimulation in memory formation and their association with amino acid tissue content and release in normal rats,” Hippocampus, vol. 22, no. 1, pp. 98–105, Jan. 2012, PMID: 20882549.