LFPsim was designed as an easy-to-use tool for modeling and computing extracellular electric field of single neuron and LFP of a population of neurons. LFPsim uses NEURON's extracellular mechanism to calculate total ionic currents from neuronal compartments. Three biophysical modeling schemas were implemented to model the extracellular activity. The field potential was calculated using PSA, LSA (Holt and Koch, 1999; Gold et al., 2006) and Resistance-Capacitance (RC) based low-pass filter techniques (Bédard et al., 2004). Extracellular potential at each time step (dt) was calculated by setting pointers to “lfp.mod” and for multiple recording points, “mea.mod.” “lfp.mod” and “mea.mod” are NMODL files used by the LFPsim to calculate extracellular potential generated by individual current sources at each time step.

Extracellular field potentials were computed from transmembrane ionic currents in the extracellular medium (Reimann et al., 2013). Transmembrane ionic currents are summed active currents generated from individual ion channels of neuronal processes and they diffuse into extracellular medium when action potential propagate in the neuron.

With NEURON, the transmembrane ionic currents in multi-compartmental models can be calculated by summing up all active currents estimated using extracellular mechanism from NEURON as implemented in LFPsim (see Equation 1). Circuit representation of intracellular and extracellular potentials is shown in Figure 1A. Transmembrane ionic currents are approximated such that the total current through a compartment is equal to the current density as the center of the compartment times the membrane area of the compartment (Segev et al., 1989; De Schutter, 2010). (1)Itransmemberane=Iionic+cm∂Vm∂t Where, I_ionic represent the ionic currents and cm∂Vm∂t represent the capacitive current.

The forward modeling schemas (Rall and Shepherd, 1968; Holt and Koch, 1999) allow modeling the extracellular field potential around neurons (Figure 1B). The computational schemas of the LFPsim involve two steps. (1) Estimation of total ionic current from detailed biophysical model of neuronal compartments. (2) Calculation of field potential at a point, P(x,y,z) from calculated ionic currents. In LFPsim, computations were estimated for each time step.

From individual current sources, the extracellular potential was computed using PSA, LSA (Holt and Koch, 1999; Gold et al., 2006) and RC low pass filter techniques, implemented within LFPsim.

Studies (Rall and Shepherd, 1968; Bédard et al., 2004; Gold et al., 2006) indicate neuropil can be modeled as an isotropic volume conductor with no capacitive effect of the medium for frequency range (0–3000 Hz). This allows the extracellular medium to be described as a purely ohmic conductivity (Plonsey, 1969; Holt, 1998).

In quasistatic approximation of Maxwell's equation, the electric field (E) and magnetic field (B) are effectively decoupled (Hämäläinen et al., 1993). ∇ ^*E = 0, the electric field in the extracellular medium can be related to extracellular potential (Φ) by E = −∇Φ. A linear relationship was assumed between the transmembrane current density (J_m) and the extracellular field potential (see Equation 2, Hämäläinen et al., 1993). σ denotes the extracellular conductivity. (2)σ∇Φ=Jm Extracellular potential generated from a neuronal compartment or segment can be generally approximated as either to a point or line source (Holt and Koch, 1999; Gold et al., 2006). Point source models assume transmembrane currents are generated from a point at the center of the neurite (see Figure 1C). Currents from each compartment across three dimensional space were computed from individual point sources. In the first schema, PSA technique was implemented in LFPsim to estimate the extracellular waveform. For a single point source, “I” denotes the transmembrane current generated from the source, “σ” denotes the conductivity of the medium and “r” denotes the distance from source to the point of measurement. Extracellular potential (Φ) at a distance “r” is computed as (3)Φ=I4πσr Considering “n” point sources in the extracellular medium, LFP was computed as, (4)ΦLFP=∑i=1n_sourcesIi4πσri Φ_LFP denotes the estimated LFP from transmembrane ionic currents.

As a second technique, LSA was implemented in LFPsim. LSA implementation assumed continuous distribution of the transmembrane currents generated from a line which passes through the axis of the neurite (see schematic in Figure 1D). LSA implementation have been known to better approximate the extracellular signal (except at very close distance, less than 1 μm away; Rosenfalck, 1969; Trayanova et al., 1990; Gold et al., 2006). LSA implementations have been previously modeled and validated on pyramidal neurons (Gold et al., 2006). Assuming line distribution of currents, extracellular potential of a line segment was estimated as (5)Φ=14πσ∫-Δs0IdsΔsr2+(h-s)2 (6)Φ=I4πσΔslogh2+r2-hl2+r2-l For “n” individual sources (line segments) in the extracellular medium, LFP was computed as in Equation (7), (7)ΦLFP=∑i=1n_sourcesIi4πσ△siloghi2+ri2-hili2+ri2-li Where, Φ_LFP denotes the calculated LFP from transmembrane ionic current, “Δs” denotes the length of single line source, “r” denotes the radial distance from the line, “h” denotes the longitudinal distance from the end of the line, l = Δs+h denotes the distance from the start of the line, “σ” denotes the conductivity of the extracellular medium (see Figure 1D for schematic) and “I” denotes the transmembrane current generated from the source.

In another schema, neuronal processes were considered as passive compartments (Figure 1E). The low pass filtering property of the extracellular medium was modeled using a simple RC (low pass) filter (Bédard et al., 2004). Extracellular potential (Φ) at a distance “d” can be calculated as (8)Φ=I e-(tEREC) For “n” sources in extracellular medium, the computation is denoted as in Equation (9), (9)ΦLFP=∑i=0n_sourcesIie-(tEREC) Where “V₀” denotes individual compartment transmembrane ionic currents, “t” denotes the time constant (t = E_rE_c). “E_C” denotes capacitance of extracellular medium, set to specific capacitance of the membrane (Johnston and Wu, 1995; Bédard et al., 2004). In this schema, homogenous capacitance was assumed throughout the extracellular space, “E_R” denoted resistance of extracellular medium and standard value was assumed (cytoplasmic resistivity of squid axon set to 0.35 Ωm; see Bédard et al., 2006, for discussion on these parameters).

The LFPsim interface includes 4 views; Morphology view window (Figure 2A) is a shape plot which helps to visualize neuron morphology and extracellular electrode location (showed in blue). Right clicking on this window allows options to the user to view the model at different angles and allows selecting LFP electrode location by selecting “LFP_electrode” menu, Voltage changes view window (Figure 2B) is a space plot to visualize voltage changes across the neuronal membrane during activity. Voltage range from −70 to +40 mV was defined as color map; computed LFP view panel (Figure 2C) allows viewing the simulated LFP trace in all three methods; simulation control (Figure 2D) includes controls for setting electrode parameters and extracellular medium properties and running simulation. LFPsim will be made publically available at ModelDB (http://modeldb.yale.edu).

Multiple Electrode Array (MEA) recordings have been extensively used in electrophysiology to study neuronal circuit function in brain slices and intact brains (BeMent et al., 1986; Mapelli and D'Angelo, 2007; Mapelli et al., 2010; Spira and Hai, 2013). Computational reconstruction of multiple electrodes was aimed to help observe, constrain and model activity at network level. In LFPsim, “MEA” like observation was modeled mathematically by spatially arranging individual electrodes in a square matrix (see Figures 1C, 3A). The tip of the individual electrode was considered as the recording point of the electrodes. The effect of saline layer interfacing between brain slice and MEA chip is not considered in the model in order to reduce the complexity (Ness et al., 2015). Field potential generated from a neuronal area was modeled using LSA (see Equations 6, 7). LFP at a point, MEA(x,y,z) was computed using Equation (10). (10)MEA(x,y,z)=ΦLFP(x,y,z)

For multi-electrode computation, a cerebellar granular layer network model was used to reproduce “MEA” LFP feature of LFPsim. The network model was simulated with center-surround excitation (Parasuram et al., 2015), with a single spike as input (in vitro behavior) through mossy fibers. The cerebellar circuit was simulated for 200 ms, input stimuli to network was provided at t = 20 ms through mossy fibers by single spike to compute in vitro post-synaptic LFP (Mapelli and D'Angelo, 2007; Diwakar et al., 2011; Parasuram et al., 2011). After loading the model in LFPsim in NEURON, multi electrode location was set by accessing “MEA” properties described in LFPsim (see Figures 3B,C). By default, an array of 4-by-4 LFP electrodes was deployed and distance between electrodes was set to 100 microns (see schematic, Figure 3A). The number of electrodes, distance and plane of electrode can be altered via the LFPsim interface.

The cerebellum granule neuron model (Diwakar et al., 2009), a detailed biophysical model consisting of 52 active cable compartments, was used to compute extracellular potential. Distribution and localization of the channels have been described previously (D'Angelo et al., 2001; Diwakar et al., 2009). The model has four excitatory (with AMPA and NMDA receptor dynamics) and four inhibitory (GABA) synapses, located at the distal end of the dendrites. A mossy fiber excites both Granule and Golgi cells of same glomeruli. In the model, the mossy fibers excited granule cells via excitatory synapses (D'Angelo et al., 2001; Nieus et al., 2006) and Golgi cell axon inhibited the granule cells via GABAergic synapse (Nieus et al., 2006).

A detailed biophysical model of the cerebellum Golgi neuron (Solinas et al., 2007) was used to construct the cerebellar network. The Golgi neuron model consisted of a soma, three dendrites and an axon. Most of the ion channels mechanisms were placed at the soma as described in Solinas et al. (2007). Synaptic model including AMPA, NMDA, GABA receptor dynamics were implemented as described in Solinas et al. (2010).

Pyramidal neuron models were adapted from Mainen and Sejnowski (1996), models of L3 Pyramidal neuron, L3 Stellate neuron and L5 Pyramidal neuron were used to study extracellular field potential computation in complex neurons. The models were stimulated by somatic current injection 70, 100, 200 pA and synaptic activations as reported in Mainen and Sejnowski (1996).

The cerebellar granular layer network model consisted of 730 Multi-compartmental granule neurons (Diwakar et al., 2009), 40 Mossy Fibers (MF) and approximately 8500 synapses to a pack 35 μm cubic slice of cerebellar cortex. The network was simulated with in vitro pattern (single spike as input at 500 Hz) to generate the N_2a and N_2b waves (Figure 5A) as reported in Mapelli and D'Angelo (2007) and in simulation by Diwakar et al. (2011) and Parasuram et al. (2011). Plasticity conditions were simulated in cerebellar granular layer by modifying intrinsic excitability of sodium channel and release probability of excitatory synapses (Nieus et al., 2006).

As an alternative model in this study, a biophysical network model of neocortex (Vierling-Claassen et al., 2010) was used, (ModelDB accession no: 141273). The model represents cortical layer II/III circuits containing three cell types: pyramidal fast spiking, regular spiking, and low-threshold spiking interneurons. The model was developed in NEURON and input to the model was varied at frequencies to generate LFP waves at 8, 32, and 56 Hz FS as mentioned in Vierling-Claassen et al. (2010).

LFP modeling schemas (PSA, LSA, and RC methods) were implemented in “extracellular_electrode.hoc.” In LFPsim, field potential calculation at each time step, dt, was coded in NMODL, “lfp.mod” was for simulating single electrode LFP simulation and “mea.mod” for multiple electrode simulation. Procedures for changing electrode recording position were implemented in “move_electrode.hoc.” “multiple_electrode.hoc” contains implementation of LFP for multiple electrode simulation. “multiple_electrode1.hoc” contains functions to set multi electrode simulation. “mea_run_then_plot.hoc,” contains the functions to plot MEA trace using a NEURON graph plot and “tool_interface.hoc” contains the GUI interface script. Details for using LFPsim for computing single extracellular potential and population LFP are included in Supplementary Material. The current version of the tool was not designed to run on parallel platforms, and therefore may not be suitable for very large scale networks.

All simulations were performed with the NEURON version 7.4 simulation environment. Single neuron and network simulations were performed on a 6-core Intel(R) Xeon(R) CPU W3670 at 3.20 GHz processor and 8GB of RAM. All code is available from ModelDB (http://modeldb.yale.edu) with accession number 190140.

Using LFPsim, extracellular field potentials from electrotonically compact cerebellar granule neuron (Diwakar et al., 2009) and morphologically complex neurons of neocortical column (Mainen and Sejnowski, 1996) were computed.

With spike inputs via 3 mossy fiber synapses, the granule neuron model generated a single action potential. Extracellular potential calculated at a point in the vicinity of the axon showed increased amplitude compared to other recording positions near the neuron (see Figure 4). The techniques, PSA and LSA approximations showed similar responses in terms of estimated extracellular signal for single neurons. The low pass filter based method exhibits reduced amplitude and width in the computed wave (see Figure 4). The nature of the single neuron extracellular potential depended on compartmental contributions and morphological details (Parasuram et al., 2011). We also observed that non-detailed models did not provide reliable field potential reconstructions (data not shown).

With the neocortical column (Mainen and Sejnowski, 1996), L3 pyramidal neuron, L3 stellate neuron, and L5 pyramidal neuron models, extracellular potential computations were similar to those in experimental studies. Somatic current injections (70–200 pA) were provided as inputs (Mainen and Sejnowski, 1996) and extracellular potential for the first spike was calculated (see Figure 4). The RC-based method showed reduced extracellular signal amplitude and width (See Figure 4) although PSA and LSA had no significant differences.

In this study, two biophysically detailed network models, implemented in NEURON were used to compute network LFPs. For structurally less complex, small scale networks, the cerebellar granular layer model was used (Parasuram et al., 2011). A network model of neocortex (Vierling-Claassen et al., 2010) was used to investigate the role of cable structures in LFP simulations on complex microcircuit models.

Three modeling methods generated extracellular field potentials of single neurons and neural microcircuits. Computed extracellular waves (based on PSA and LSA) showed linear drop in amplitude when the recording electrode moved away from the neuronal process (see Figures 7A,E). LFP wave computed from (PSA and LSA) also showed the similar amplitude and shape for recording sites farther than 120 microns from the neurite (see Figure 7B).

Extracellular potentials of a granule neuron, computed when the electrode was in the vicinity to the soma (Figure 7D, LSA in red, PSA in green, Simple RC in blue) showed similar behavior. LSA and RC based LFP waveforms showed similar amplitude but the PSA generated LFP wave showed slightly reduced amplitude compared to the other techniques.

Attenuation properties of the extracellular medium was studied on a multi compartmental model of the granule neuron (Diwakar et al., 2009) by varying extracellular resistivity from 0.24 to 0.42 Ωm (Nicholson and Freeman, 1975; Goto et al., 2010). Extracellular wave amplitude showed a linear decrease, when the resistivity of extracellular medium increased (see Figure 7C).

LFP algorithmic implementations were analyzed using RAM model of order analysis (Hartmanis, 1971). PSA-based estimations of extracellular electric potential of a single current source at a given point for a single time step (dt) used in 21 unit operations (each arithmetic operation was counted as one unit in order analysis methods), LSA used 70 units and RC filter used 21 units. For a neuron model with “n” compartments, 70^*n units for LSA to calculate the electric potential for single time step. For large “n” on network models, the algorithms will run in the order of TS, O (TS) where “TS,” denotes the total time steps in the simulation.

Computation time and memory requirements to simulate large scale network models of LFPsim was analyzed. Real time measurements for computed LFPs were plotted corresponding to a 100 ms of simulation (See Figure 7F). Augmenting granule neuron population sizes indicated a linear increment in computing time when network model was scaled by increasing number of current sources (Figure 7F). The amount of RAM used to compute LFP for a model did not vary much throughout the simulation (Figure 7G).

Using multiple spaced virtual recording electrodes, in vitro “MEA” recording of granular layer and spatial attenuation of cerebellar N_2a and N_2b LFP waves (see Figure 3D) was simulated in a 400 × 400 μm area of granule neurons with 100 μm spacing among electrodes. LFP computed from peripheral electrodes showed decreased amplitude and width of N_2a and N_2b waves (Diwakar et al., 2011) in comparison to the center of MEA (Figure 3D). Array-based representation also reproduced the inherent amplitude decay modeled as a characteristic of the extracellular space.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.